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wavs.  Either  a  complex  higher  order  solution  scheme  is  employed  or  a  fine 
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spacing  is  constructed  over  the  entire  solution  domain  so  that  low-order 
schemes  can  be  used.  From  a  programming  and  computational  efficiency  view¬ 
point,  low-order  schemes  are  more  desirable.  However,  increasing  the  number 
of  computational  points  is  often  an  inefficient  approach  since  portions  of 
the  domain  may  contain  small  gradients  of  the  solution  variables  and  a  fine 
grid  is  not  required  there. 

An  approach  that  can  sometimes  be  taken  which  reduces  numerical  error 
while  still  allowing  for  low-order  schemes  is  the  use  of  an  adaptive  grid. 

Such  a  grid  moves  with  the  developing  solution  so  that  a  small  grid  spacing 
occurs  in  regions  where  the  solution  gradients  are  large  with  a  large  spacing 
occurring  in  regions  where  the  solution  is  more  uniform.  Grid  movement  al¬ 
gorithms  can  be  generated  either  by  employing  a  variational  principle  or  by 
allowing  the  points  to  attract  and/or  repel  each  other  in  a  manner  similar  to 
the  way  electrical  charges  behave.  In  either  case,  some  measure  of  the  error 
is  equidistributed  over  the  grid.  Although  an  adaptive  grid  normally  consists 
of  a  fixed  number  of  points,  variable  node  grids  can  also  be  constructed  by 
adding  and/or  deleting  points  in  particular  regions  so  as  to  equidistribute 
the  numerical  error. 

•This  investigation  has  centered  around  the  assessment  of/ the  usefulness 
of  adaptive  grids  in  the- numerical  solution  of  free  surface  hydrodynamics  with 
particular  problems  expected  to  occur  in  the  use  of  adaptive  grids  in  hydro- 
dynamic  mode lingy being  noted,  e.g.,  the  need  for  interpolation  of  the  bottom 
topography  on  a  moving  grid.  Methods  of  grid  generation,  as  well  as  the  im¬ 
plementation  of  adaptive  grids,  have  been- considered  for  both  -*he‘  finite  dif¬ 
ference  and  the  finite  element  solution  methods. 

In  general,  it  can  be  concluded  that  adaptive  grids  offer  significant 
potential  for  more  accurate  solutions  at  less  cost  when  modeling  the  behavior 
of  surges,  hydraulic  jumps,  and  concentration  fronts,  i.e.,  perturbation 
problems  containing  only  a  few  high  gradient  regions  in  the  computational  do¬ 
main  at  any  time.  However,  for  problems  such  as  the  propagation  of  short  waves 
in  coastal  regions,  where  many  wavelengths  occur  in  the  physical  domain  of  in¬ 
terest  at  a  particular  time,  adaptive  grids  offer  little  or  no  advantage  since 
gradients  occur  over  essentially  the  complete  domain. 

Some  of  the  greatest  potential  for  adaptive  grid  techniques  in  hydro- 
dynamic  modeling  lies  in' the- use  of  -such  techniques  to  numerically  generate 
fixed  grids  in  tidal  circulation  studies.  Such  studies  normally  involve  estu¬ 
aries  and  coastal  areas  containing  navigation  channels  requiring  fine  grid 
resolution.  Adaptive  grid  techniques  employing  adaption  to  the  water  depth 
should  be  of  great  use  in  developing  grids  for  such  problems,  and  research  ef¬ 
forts  in  this  direction  are  strongly  encouraged. 


_ Unclassified _ 

SECURITY  CL  A  SSI  Etc  ATlON  OF  THIS  P  AGEfWTien  Omtm  Entered) 


PREFACE 


The  study  reported  herein  was  conducted  during  the  period  February  1984 
to  '<_•,) t<.  Tiber  1984  by  the  Hydraulics  1  iborutory  of  the  US  Army  Engineer  Water¬ 
ways  Experiment  Station  (WES)  under  the  general  supervision  of  Messrs.  H.  B. 
Simmons,  former  Chief  of  the  Hydraulics  Laboratory,  and  F.  A.  Herrmann,  .Jr., 
Chief  of  the  Hydraulics  Laboratory,  and  M.  B .  Boyd,  Chief  of  the  Hydraulic 
Analysis  Division  (HAD).  The  study  was  funded  by  Department  of  the  Army  Proj¬ 
ect  4A0iS  1 10 1 A9 ID,  "In-House  Laboratory  Independent  Research,"  sponsored  by  the 
Assistant  Secretary  of  the  Army. 

Dr.  B.  H.  .Johnson,  HAD,  directed  the  study  and  the  preparation  of  the 
report,  with  assistance  from  the  coauthors  Dr.  J.  F.  Thompson  of  the  Depart¬ 
ment  of  Aerospace  Engineering  at  Mississippi  State  University  and  Dr.  A.  J. 
Baker  of  the  Department  of  Engineering  Science  and  Mechanics  at  the  University 
of  Tennessee. 

Commanders  and  Directors  of  WES  during  the  conduct  of  this  study  and  the 
preparation  and  publication  of  this  report  were  COL  Tilford  C.  Creel,  CE,  and 
COL  Robert  0.  Lee,  CE.  Technical  Director  was  Mr.  F.  R.  Brown. 


CONTENTS 


Page 


PREFACE  .  1 

CONVERSION  FACTORS,  US  CUSTOMARY  TO  METRIC  (SI)  UNITS 

OF  MEASUREMENT  .  3 

PART  I:  INTRODUCTION .  4 

Types  of  Numerical/ Free  Surface  Hydrodynamic  Models  .  5 

Solution  Techniques  .  .  .  6 

Adaptive  Grids . II 

PART  II:  IMPLEMENTATION  OF  AN  ADAPTIVE  GRID .  15 

Finite  Difference  Modeling  ....  .  15 

Finite  Element  Modeling  .  ..  .  17 

PART  III:  GENERATION  OF  ADAPTIVE  GRIDS .  25 

Fixed  Number  of  Grid  Points  . .  25 

Variable  Number  of  Grid  Points .  40 

PART  IV:  APPLICABILITY  OF  ADAPTIVE  GRIDS  IN  HYDRODYNAMIC  MODELING  .  .  43 

Basic  Concerns  . . 43 

Hydrodynamic  Applications  .  47 

PART  V:  SUMMARY .  55 

Conclusions .  55 

Recommendations . 56 

REFERENCES .  58 

FIGURES  1-15 

APPENDIX  A:  HYDRODYNAMIC  EQUATIONS  AND  APPROXIMATIONS  .  A1 


CONVEKSION  FACTORS,  US  CUSTOMARY  TO  METRIC  (SI) 
UNITS  OF  MEASUREMENT 


A  DISCUSSION  OF  ADAPTIVE  GRIDS  AND  THEIR  APPLICABILITY 


IN  NUMERICAL  HYDRODYNAMIC  MODELING 


PART  I:  INTRODUCTION 


1.  When  attempting  to  determine  the  hydrodynamics  of  water  bodies  such 
as  rivers,  estuaries,  etc.,  the  US  Army  Corps  of  Engineers  uses  one  or  perhaps 
a  combination  of  three  approaches — field  investigations,  physical  models,  and 
mathematical  or  numerical  models.  FieLd  investigations  may  reveal  what  pres¬ 
ently  occurs  in  a  water  system  but  cannot  predict  what  will  result  from 
changes  due  to  new  inputs  to  the  system.  Depending  upon  their  complexity, 
physical  models  can  require  significant  investments  of  capital  and  long  con¬ 
struction  and  test  periods;  whereas,  depending  upon  approximations  made  to  the 
governing  equations  of  motion  and  the  solution  technique  employed,  numerical 
models  can  often  provide  relatively  low  cost  and  highly  flexible  models  for 
some  problems. 

2.  Historically,  the  majority  of  numerical  hydrodynamic  models  have  em¬ 
ployed  the  method  of  finite  differences  to  solve  the  equations  of  motion  on  a 
rectangular  grid  with  a  uniform  and  fixed  grid  spacing.  As  a  result,  irregular 
boundaries  can  only  be  represented  in  a  "stair-stepped"  fashion.  In  addition, 
with  a  fixed  and  uniform  grid  spacing,  numerical  errors  normally  occur  in 
regions  where  large  solution  gradients  develop  unless  the  grid  spacing  is  small 
enough  to  accurately  resolve  the  developing  gradients.  Such  grids  often  con¬ 
tain  so  many  computational  points  as  to  make  computations  uneconomical.  Over 
the  past  several  years,  numerous  efforts  have  been  made  to  address  these  prob¬ 
lems.  The  use  of  boundary-fitted  grids  in  finite  difference  models,  as  well 

as  the  development  of  finite  element  hydrodynamic  models,  are  examples  of  ap¬ 
proaches  for  handling  irregular  geometry,  as  well  as  grid  resolution  problems. 
If  the  locations  of  solutions  gradients  are  known  a  priori.  However,  the  loca¬ 
tions  of  these  gradients  are  not  normally  known  beforehand,  especially  when 
modeling  moving  disturbanc les ,  e.g.,  surges  resulting  from  hydropower  peaking 
operations . 

3.  For  the  past  two  years,  a  significant  amount  of  research  has  been 
ongoing  in  the  aerodynamics  community  on  the  development  of  adaptive  or 
moving  grid  techniques  to  address  grid  resolution  problems.  This  report  is 


4 


an  attempt  to  assess  the  usefulness  of  such  techniques  in  the  numerical  model¬ 
ing  of  free  surface  hydrodynamic  problems. 


Types  of  Numer ica 1/Free-Sur face  Hydrodynamic  Models 

4.  Numerical  hydrodynamic  models  can  differ  widely,  depending  upon 
such  things  as  the  solution  technique  applied  to  the  governing  differential 
equations  representing  the  physical  processes,  the  assumptions  made  in  the 
derivation  of  the  governing  equations,  whether  the  phenomena  are  steady  or 
time-varying,  and  the  spatial  dimensionality  considered,  with  perhaps  the  spa¬ 
tial  dimensionality  being  the  most  commonly  used  delineator. 

5.  If  the  complete  three-dimensional  (3-D)  equations  of  motion  are  in¬ 
tegrated  over  a  cross  section,  one-dimensional  ( 1— D)  models  result.  Such 
models  are  commonly  applied  in  computing  river  hydrodynamics,  e.g. ,  the  compu¬ 
tation  of  floods.  Averaging  over  either  the  depth  or  the  width  results  in 
two-dimensional  (2-D)  models.  Vertically  averaged  models  are  applicable  for 
the  computation  of  nearly  horizontal  flow  in  relatively  shallow  and  well-mixed 
bodies  of  water,  whereas  laterally  averaged  models  are  appropriate  when  deal¬ 
ing  with  relatively  narrow  and  deep  bodies  of  water  experiencing  vertical 
stratification  of  the  water  density. 

ft.  Even  though  a  free  surface  exists  on  open  bodies  of  water,  some 
modelers  have  treated  the  surface  as  a  rigid  lid  when  very  little  motion  of 
the  free  surface  occurs.  The  surface  then  becomes  in  essence  a  solid  boundary 
and  the  normal  component  of  the  velocity  must  be  zero.  In  addition,  the  pres¬ 
sure  can  no  longer  be  prescribed  at  the  surface  but  rather  must  be  computed. 
The  pressure  boundary  condition  then  takes  the  form  of  a  derivative  boundary 
condition,  i.e.,  a  Neumman  condition  as  opposed  to  a  Dirichlet  condition  in 
the  true  free-surface  case.  All  hydrodynamic  modeling  discussed  in  this  re¬ 
port  treats  the  surface  as  being  free  to  move  so  that  free-surface  waves,  e.g. 
tidal  waves  in  estuaries,  are  free  to  be  computed.  In  addition,  only  1-  and 
2-D  problems  are  considered  since  little  research  on  3-D  adaptive  grids  has 
been  conducted  to  date. 

7.  For  completeness,  a  discussion  of  common  assumptions  made  in  free- 
surface  hydrodynamic  models,  e.s?.,  a  hydrostatic  pressure,  as  well  as  the 
governing  equations  for  both  1-  and  2-D  models,  are  presented  in  Appendix  A. 
Various  types  of  problems  to  which  such  models  are  applied  and  suggestions  on 


the  use  of  adaptive  grid  techniques  in  the  application  of  those  models  are  the 
major  subjects  to  be  discussed  in  the  text. 

Solution  Techniques 

8.  The  governing  hydrodynamic  equations  presented  in  Appendix  A  are  non¬ 
linear  partial  differential  equations,  which  in  a  strict  mathematical  sense 
are  classified  as  being  of  the  parabolic  type.  However,  outside  the  boundary 
layer  the  equations  exhibit  a  strong  hyperbolic  or  wave  character  due  to  the 
dominance  of  the  convective  terms  and  thus  are  often  considered  as  being  of 

the  hyperbolic  type.  In  any  case,  because  of  the  nonlinearity,  analytical 
solutions  do  not  generally  exist  and  one  must  resort  to  numerical  methods  to 
obtain  an  approximation  of  the  continuous  solution  of  the  differential  equa¬ 
tions.  Such  methods  consist  primarily  of  the  use  of  either  finite  differences 
or  finite  elements. 

Finite  element  method 

9.  In  the  finite  element  approach,  the  field  is  divided  into  smaller 
regions  of  convenient  shapes,  such  as  triangles  or  quadrilaterals,  and  the 
solution  is  approximated  on  each  element  by  interpolation  from  nodal  values  on 
the  element.  Using  a  variational  principle,  or  a  weighted-residual  method, 
the  governing  partial  differential  equations  are  then  transformed  into  finite 
element  equations  governing  each  isolated  element.  These  local  equations  are 
then  collected  to  form  a  global  system  of  ordinary  differential  (in  time),  or 
algebraic,  equations  including  a  proper  accounting  of  boundary  conditions. 

The  nodal  values  of  the  dependent  variables  are  then  determined  from  solution 
of  this  matrix  equation  system. 

10.  Despite  its  ability  to  resolve  complex  geometry  in  the  computational 
domain,  the  finite  element  method  suffers  from  the  fact  that  it  is  relatively 
more  complicated  and  expensive  to  program.  In  addition,  most  existing  hydro- 
dynamic  finite  element  models  appear  to  require  significantly  more  computational 
effort  per  time-step.  However,  this  may  be  attributed  to  the  techniques  com¬ 
monly  employed  by  finite  element  modelers  in  solving  the  matrix  equation  system 
rather  than  the  method  itself. 

Finite  difference  method 

IL.  The  majority  of  numerical  hydrodynamic  models,  whether  they  be  1-, 

2-,  or  3-D,  employ  the  use  of  finite  differences  to  obtain  solutions  of  the 


riverains  equations  of  fluid  motion.  In  the  finite  difference  method,  the  do¬ 
main  of  tiie  independent  variahLes  is  replaced  by  t  finite  set  of  points  re¬ 
ferred  to  as  net  or  mesh  points.  One  tiien  seeks  to  determine  approximate 
alue^  for  the  desired  solutions  at  these  points.  The  values  at  the  mesh 
points  are  require-  to  satisfy  difference  equations  that  are  usually  obtained 
by  replacing  partial  derivatives  by  partial  difference  quotients.  The  result¬ 
ing  set  of  simultaneous  algebraic  equations  is  then  solved  for  the  values  of 
the  solution  at  the  mesh  points.  If  the  boundaries  do  not  coincide  with  mesh 
points,  the  finite  difference  approach  applied  to  the  equations  in  a  Cartesian 
coordinate  system  requires  interpolation  between  mesh  points  to  represent 
boundary  conditions.  However,  through  boundary-fitted  coordinate  transforma¬ 
tions  irregular  boundaries  can  he  accurately  handled  while  stilL  making  use  of 
the  simplicity  of  finite  differences  to  obtain  solutions.  The  most  general  of 
such  transformations,  which  is  discussed  in  more  detail  in  the  next  section, 
is  a  method  developed  by  Thompson  et  al .  (1974),  which  generates  curvilinear 
coordinates  as  the  solution  of  two  eLliptic  partial  differential  equations. 
Crids  employed 

12.  Whether  the  finite  element  or  finite  difference  method  is  employed, 
the  Dody  of  water  being  modeLed  must  be  covered  with  a  discrete  set  of  points. 
These  poinis  constitute  the  numerical  grid  upon  which  a  solution  of  the  govern¬ 
ing  equations  is  sought,  and  depending  upon  the  placement  of  these  points  and 
the  relative  magnitude  of  solution  gradients,  numerical  errors  invariably  occur 
As  previou.,^/  noted,  finite  element  grids  generally  consist  of  a  mixture  of 
triangles  and  quadrilaterals.  Within  the  framework  of  finite  difference  solu¬ 
tions,  there  have  been  several  attempts  to  obtain  solutions  on  curvilinear,  as 
welL  as  stretched  rectangular,  grids  in  order  to  better  handle  irregular  bound¬ 
aries,  as  weLI  as  to  improve  grid  resolution  problems  and  thus  to  decrease  the 
numerical  error. 


13.  Wanstrath  (1977)  developed  a  vertically  averaged  numerical  model 
for  computing  storm  surges  at  coastlines  through  the  use  of  conformal  mapping. 
A  spatial  region  of  prototype  space  is  conformal ly  mapped  into  a  rectangle  in 
a  mathematical  image  plane.  To  provide  an  evenly  spaced  computing  grid,  an 
independent  stretching  of  the  coordinates  is  then  performed. 

14.  The  Waldrop  and  Ttfom  (’1976)  hydrodynamic  model  was  developed  to 

analyze  hr  t  water  discharges  from  power  plants  er*  uses  an  orthogonal  curvi¬ 
linear  griil  that  is  generated  i !  gehr  i  i  ■  i  i  ;  v  ,  !  y  interpolation.  As  with 


the  grids  generated  by  Wanstrath  through  conformal  mapping,  these  grids  lack 
generality  in  that  boundary  point  selection  and  grid  spacing  are  not  arbitrary. 
In  addition,  interior  bodies  are  not  allowed  in  either  the  Wanstrath  or  the 
Waldrop  models.  More  recently,  Blumberg  and  Herring  (1983)  have  developed  a 
3-1)  hydrodynamic  model  which  also  handles  the  horizontal  dimensions  through 
the  use  of  an  orthogonal  curvilinear  grid. 

15.  In  addition  to  Wanstrath,  others  such  as  Waldrop  and  Farmer  (1976) 
and  Butler  (1978)  have  employed  the  use  of  independent  stretching  of  the  hori¬ 
zontal  coordinates  to  improve  grid  spacing.  For  example,  the  stretching  trans¬ 
formation  in  Waldrop's  buoyant  plume  computations  is 


X 


1  - 1  x 

—  tan 

C1  1 2 


Y 


1  - 1  y 

-  tan 
cn  c . 


(1) 


thus  if  Cj  =  ;/2  ,  then  at  x  =  •>>  ,  X  =  1.9  .  Even  increments  of  AX  thus 
produce  a  close  spacing  of  grid  points  near  the  plume  entrance,  where  the  best 
resolution  is  desired,  and  yet  extend  the  region  of  computation  far  from  the 
origin  so  that  boundary  conditions  can  be  more  easily  specified.  Other  re¬ 
searchers  have  employed  different  stretching  functions,  e.g.,  Butler  employs 
an  exponential  form  of  stretching. 

16.  Vertically  averaged  hydrodynamic  models  developed  by  Johnson  (1980) 
and  Siieng  and  flirsli  (1984)  are  examples  of  models  employing  a  general  nonor- 
thogonal  curvilinear  grid.  Through  coordinate  transformations,  irregular 
boundaries  and  variable  grid  spacing  can  be  more  accurately  handled  while  still 
making  use  of  the  simplicity  of  finite  differences  to  obtain  solutions.  The 
extremely  general  coordinate  transformation  employed  results  from  a  method  de¬ 
veloped  by  Thompson  et  al.  (1974),  which  generates  curvilinear  coordinates  as 
the  solution  of  two  elliptic  partial  differential  equations  with  Diriehlet 
boundary  conditions,  one  coordinate  being  specified  as  constant  on  the  bound¬ 
aries,  and  a  distribution  of  the  other  specified  along  the  boundaries.  No 
restrictions  are  placed  on  the  irregularity  of  the  boundaries,  and  fields  con¬ 
taining  multiple  bodies  or  branches  can  be  handled  as  easily  as  simple  geom¬ 
etries.  ilega rd 1  ess  of  the  shape  and  number  of  bodies  and  regardless  of  the 
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>  pae  i  in.’,  >!  ci'orJin;iU‘  lines,  all  numerical  computations,  both  to  generate  the 
coordinate  system  and  to  subsequent  1 y  solve  the  fluid  flow  equations,  are  done 
on  a  rectangular  grid  with  square  mesh. 

! 7 .  Since  t  lie  boundary-lilted  coordinate  system  has  a  coordinate  line 
coincident  with  all  boundaries,  all  boundary  conditions  can  be  expressed  at 
grid  points;  and  normal  derivatives  can  be  represented  using  only  finite  dif¬ 
ferences  between  grid  points  on  coordinate  lines.  Mo  interpolation  is  needed, 
even  though  the  coordinate  system  is  not  orthogonal  at  the  boundary. 

18.  As  will  be  discussed  later  in  connection  with  generation  techniques 
for  adaptive  grids,  a  logical  choice  of  the  elliptic  generating  system  is 
Poisson's  equation.  Thus  for  a  domain  such  as  illustrated  in  Figure  1,  the 
basic  problem  is  to  solve 


+ 

XX  vv 


(2) 


n  +  n  =  Q 

xx  yy 


with  boundary  conditions 


=  y*. 

V)  , 

•;  =  const. 

on 
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:  i  ’ 
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=  const, 

on 

r3 

=  yx. 

y)  , 

a  =  const. 
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5  ; 

II 

X 

•  < 

=  const  , 

on 

1  7 

j  '.y. 

y )  i 

=  const. 
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f  =  const, 
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=  Vx’ 

v)  , 

=  const. 

on 
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6  ’ 

■  =  ng(x,v)  , 

A  =  const. 

on 
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The  functions  P  and  Q  may  be  chosen  to  cause  the  coordinate  lines  to  con¬ 
cent  rati-  as  desired.  One  form  of  these  functions  incorporated  by  Thompson  is 
that  of  decaying  exponentials. 

19.  Since  all  numerical  computations  are  to  be  performed  in  the  rec¬ 
tangular  transformed  plane,  it  is  necessary  to  interchange  the  dependent  and 
independent  variables  in  Fqu.it  ion  P  .  To  ace. v -pi  isl;  the  transformation,  the 
following  expressions  are  utilized 


The  system  defined  by  Equation  2  then  becomes 


ux  -  2dx„  -  yx  +  J  (Px  +  Ox  )  =  0 
>..n  tin  t.  n 


(5) 


'ty„,  -  2Ey  -  yy  +  J  (Py  +  Qy  )  =  0 
>,t,  fr)  •'nn  n 


where 

2  2 

a  =  x  +  y 
n 

B  =  x,xr|  +  y.y^  (6) 

2  2 

Y  =  x.  +  vr 

J  =  Jacobian  of  the  transformation  =  x,y  -  x  y 

with  transformed  boundary  conditions  from  Equation  3. 

20.  Although  the  new  system  of  equations  is  more  complex  than  the  origi 
nal  system,  the  boundary  conditions  are  specified  on  straight  boundaries  and 
the  coordinate  spacing  in  the  transformed  plane  is  uniform.  Computationally, 
these  advantages  outweigh  disadvantages  resulting  from  the  extra  complexity 

of  the  equations  to  be  solved. 

21.  The  rectangular  transformed  grid  is  set  up  to  be  the  size  desired 
for  a  particular  problem.  Since  the  actual  values  of  F  and  n  are  irrele¬ 
vant  in  the  transformed  plane,  the  n-lines  are  assumed  to  run  from  1  to  the 
number  of  n-lines  desired  in  the  physical  plane.  Likewise,  the  i-lines  are 
numbered  1  to  the  number  of  -lines.  The  grid  spacing  in  both  the  r  and  n 
directions  of  the  transformed  plane  is  taken  as  unity  since  these  spacings 
cancel  out  of  all  difference  expressions. 

22.  The  same  procedure  can  be  extended  to  regions  that  are  more  than 
doubly  connected,  i.e.,  have  more  than  two  closed  boundaries,  or  equivalently. 


tii. 111  Mir  body  witliin  ,1  single  miter  boundary.  A  river  reach  containing 
more  than  one  island  is  an  example. 

2  5.  As  a  i  inal  note,  both  conformal  and  orthogonal  coordinate  systems 
~>pi  •  i  i:  isrs  ot  boundary-fitted  coordinates  as  generated  from  elliptic 

systems.  As  will  be  seen  later,  the  most  promising  approach  for  generating 
adaptive  grids,  which  is  based  upon  variational  principles,  is  a  natural  ex¬ 
tension  ot  Thompson's  work  discussed  above. 

Adaptive  Grids 

Basic  idea 

24.  Tn  its  simplest  form,  an  adaptive  grid  is  considered  to  be  a 
boundary-conforming  curvilinear  coordinate  system  which  is  numerically  gen¬ 
erated  such  that  the  location  of  the  organized  set  of  points  formed  by  the 
intersections  of  the  coordinate  Lines  is  dictated  by  the  physics  of  the 
problem  being  solved.  Ultimately  then,  the  grid  points  must  be  congregated 
in  regions  of  high  activity,  i.e.  ,  Large  solution  gradients,  and  dispersed 

in  those  regions  of  the  physical  domain  where  smalL  gradients  of  the  solution 
exist.  A  good  example  of  a  1-D  adaptive  grid  used  in  the  solution  of  the 
nonlinear  Burger's  equation  is  presented  in  Figure  2.  As  shown  in  the  x-t 
plot  of  grid  points,  as  the  moving  disturbance  generated  by  the  nonlinearity 
becomes  stronger,  grid  points  are  concentrated  more  and  more  in  the  high 
gradient  region.  In  other  words,  the  physics  of  tiie  solution  being  generated 
dictates  where  the  grid  points  are  located.  A  good  discussion  of  adaptive 
grids  can  be  found  in  Chapter  XI  of  lhompson,  Warsi,  and  Mast  in  (1984),  and  a 
recent  survey  has  been  given  by  Thompson  (1984). 

Need  for  adaptive  grids 

25.  When  using  numerical  techniques  to  compute  solutions  of  partial 
differential  equations,  continuous  problems  are  cast  in  a  discretized  form. 

As  a  result,  numerical  errors  are  generally  always  present.  A  simple  illustra 
tion  of  this  is  provided  by  the  commonly  used  implicit  forward  in  time  and 
centered  in  space  finite  difference  (FI'CS)  scheme  applied  to  the  transport 
equat  ion 


where 


<p  =  any  concentration 
u  =  fluid  velocity 

x,t  =  spatial  distance  and  time,  respectively 
Applying  the  implicit  FTCS  method  yields 


,n+l  *,n 

*_i _ V;_i 

At 


2Ax 


n+1 

i+1 


( :>u) 


n+1 

i-1 


] 


(8 


Evaluating  each  term  by  expanding  in  a  Taylor  series  about  (x,t)  and  substi- 

2  3 

tuting  into  the  difference  equation  yields,  to  0(At  ,Ax  ),  the  differential 
equation 


2 

3<J>  _  S^u  9  f> 

3t  "  "  9x~  +  “e  ~2 
3x 


(9 


whe  re 


a 

e 


2  2 

A_tu_  _  A x_  3u 

2  2  ax 


(10 


Since  Equation  7  is  of  hyperbolic  form,  for  accuracy  the  numerical  domain  of 
influence  should  not  exceed  that  imposed  by  the  differential  equation,  i.e.. 

At  u  <  Ax  (11 


Using  this,  the  expression  for  can  be  written  as 


a 

e 


Ax  3u  \ 

u  3x  J 


(12 


Since  a  represents  a  diffusion  coefficient,  physically  its  value  can  never 
be  less  than  zero.  Therefore 


i.e.,  1  low  gradients  must  be  sufficiently  resolved,  which  in  turn  is  accom¬ 
plished  b\  reducing  tlie  grid  spacing  Ax  .  Obviously,  one  way  to  accomplish 
this  is  to  increase  the  number  of  grid  points.  However,  if  significant  por¬ 
tion;  of  i  lie  physical  domain  contain  relatively  small  solution  gradients  at 
any  particular  time,  much  of  the  computational  effort  is  wasted. 

26.  in  addition  to  numerical  errors  associated  with  stability  problems, 
implicit  schemes  that  are  unconditionally  stable  still  generally  contain 
numerical  errors  as  a  result  of  the  grid  spacing  that  influence  the  accuracy 
of  the  solution.  As  an  example,  consider  the  simplified  form  below  of  the 
long-period  water  wave  equations. 


cu 


+  h  :  -  =  o 

><t  >x 


04) 


Y  +  g  'P-  =  0 

•  it  ,<X 


where 

li  =  total  water  depth 

!  =  surface  elevation  relative  to  some  datum 
If  the  difference  equations  are  written  as 


n+1 
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!■  ‘  '■ 


n+l  n+1 
i+1  "  Ui-1 


=  0 


(15) 


n+!  n  , 
i  .  -  u  .  + 


't  (  n+I 


n+1 


Ax  V  ' i+1 


v  i-l 


=  0 


it  can  he  shown  that  the  computations  are  unconditionally  stable.  However,  as 
illustrated  in  Figure  5,  the  accuracy  of  both  the  amplitude  and  the  phase  of 
the  computed  wave  is  highly  dependent  on  the  number  of  grid  points  within  a 
wavelength  so  that  increasing  Lite  number  of  grid  points  obviously  will  in¬ 
crease  the  accuracy  of  the  computations. 

27.  Two  alternatives  to  increasing  the  number  of  grid  points  to  allevi¬ 
ate  stability  and/or  accuracy  problems  exist.  Higher  order  solution  schemes 
could  be  constructed.  However,  such  scheme:-  are  invariably  more  complex 
and  difficult  to  program,  especially  near  bound  tries.  In  addition,  unless 


filtering  schemes,  such  as  the  flux-corrected  transport  method  of  Zalesak 
(1979),  are  employed  higher  order  schemes  can  often  accentuate  the  solution 
"wiggles"  resulting  from  the  cell  Reynolds  number  being  greater  than  2.  In¬ 
stead  of  improving  the  solution,  a  deterioration  may  occur  even  though  the 
solution  is  more  accurate  from  an  order  of  truncation  error  analysis.  The 
second  alternative  is  to  continue  using  simple  lower  order  solution  schemes 
but  to  allow  the  grid  points  to  move  so  that  the  grid  spacing  becomes  smaller 
in  regions  where  solution  gradients  are  developing,  thereby  increasing  the 
resolution  where  needed  and  thus  decreasing  numerical  errors.  In  other  words 
compute  the  grid  as  part  of  the  overall  solution.  The  method  of  character¬ 
istics,  as  applied  in  the  solution  of  hyperbolic  equations,  can  be  considered 
as  a  type  of  adaptive  grid. 


PARI’  r  I  :  IMPLEMENTATION  OF  AN  ADAPTIVE  OKU) 


28.  When  employing  a  moving  or  adaptive  grid  to  obtain  more  accurate 

:  umericai  solutions  of  partial  differential  equations  with  a  limited  number  of 
points,  a  major  question  to  be  resolved  is  how  the  grid  is  to  be  generated 
as  time  marches  on.  Various  grid  generation  techniques  will  be  discussed 
later.  Irrespective  of  how  the  new  grid  is  generated,  a  separate  question 
concerns  the  actual  implementation  of  the  time-dependent  grid  within  the  over¬ 
all  solution  scheme  for  solving  the  differential  equations  governing  the 
physical  problem. 

F inite  _l)if  f erence  Model  ing 

29.  The  use  of  adaptive  grids  in  the  solution  of  the  hydrodynamic  equa¬ 
tions  of  motion  presented  in  Appendix  A  by  the  finite  difference.  (FD)  method 
is  best  based  upon  the  use  of  boundary-fitted  curvilinear  coordinates.  As 
discussed  in  paragraph  19,  since  all  numerical  computations  to  solve  the  hy¬ 
drodynamic  equations  are  to  be  performed  in  a  rectangular  transformed  plane 
(Figure  1),  it  is  necessary  to  interchange  the  independent  variables  (x,y)  in 
the  physical  space  and  the  (f,  ,ri)  variables  in  t he  transformed  space.  Equa¬ 
tion  4  is  used  to  accomplish  such  a  transformation .  The  resulting  transformed 
equations  have  the  curvilinear  coordinates,  r  and  n  ,  as  independent  variables, 
so  all  computations  can  be  done  on  the  rectangular  transformed  region  regard¬ 
less  of  the  shape  of  the  physical  region.  It  might  be  noted  that  in  addition 

to  transforming  the  independent  variables  the  components  of  the  velocity  vector 
can  also  be  transformed  to  yield  equations  involving  either  the  contravariant 
or  covariant  velocities.  This  approach  has  been  taken  by  Sheng  and  Hirsh 
(1984)  in  their  development  of  a  vertically  averaged  flow  model. 

30.  A  simplified  form  of  the  vertically  averaged  hydrodynamic  equations 
given  in  Appendix  A,  where  the  equations  have  been  transformed  by  using  Equa¬ 
tion  4,  is  given  below 

Continuity:  ~  ^  y  (hu) .  -  v.(lm)  +  x_(hv)  -  x  (hv)  .  =  0  (16) 

)t  Jin  n  n  n 
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x-Momentum : 


—  +  ,  (uy  -  vx  )u.  +  ,  (vx.  -  uy. )n 
it  .1  n  n  J  ■.  n 


+  ,  yng*r 


J  y>g*n 


2  2  1  /“ 
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y-Momentum:  -V  +  ^  (uy  -  vx  )v  +  ^  (vxr  -  uy,)v 
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+  J  xfg*n  "  J  Xr^r 


2  2 
m(u  +  v  ) 


All  symbols  are  defined  in  Appendix  A. 

31.  Assuming  that  the  total  number  of  grid  points  remains  the  same,  an 
adaptive  grid  could  be  implemented  with  a  solution  of  the  equations  as  written 
above.  The  equations  would  first  be  solved  on  the  initial  grid.  For  example, 
the  initial  grid  might  have  been  generated  using  Thompson's  coordinate  genera¬ 
tion  code  (Thompson  1983a).  Information  from  this  solution  would  then  be  used 
to  generate  a  new  grid  satisfying  certain  criteria  to  be  discussed  later.  In 
order  to  continue  the  time  marching  of  the  solution  on  the  new  grid,  initial 
values  of  the  solution  variables  are  required  at  each  grid  point.  These  would 
be  obtained  through  an  interpolation  of  the  computed  solution  at  the  current 
time  level  from  the  old  grid.  Depending  upon  the  accuracy  of  the  interpola¬ 
tion  used,  errors  can  be  introduced.  In  any  event,  the  interpolation  will  re¬ 
sult  in  increased  computer  time.  In  addition,  with  this  implementation,  the 
generation  of  the  grid  will  always  lag  the  flow  solution.  Obviously,  it  is 
not  necessary  to  recompute  the  grid  each  time-step  if  solution  gradients  aren't 
developing  rapidly. 

32.  A  more  natural  implementation  of  an  adaptive  grid  in  a  finite  dif¬ 
ference  solution  involves  transforming  not  only  the  spatial  (x,y)  derivatives 
but  also  the  time  derivative.  With  this  approach,  no  interpolation  is  required 
As  in  the  derivation  of  the  transformation  relations  given  by  Equation  4,  it 
can  be  shown  that  the  time  derivative  in  the  (x,y)  system  can  be  transformed 

as  below 
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<  h  t  arbitrary  variable.  Note  that  while  the  time  derivative  on  the 

lei t -hand  side  is  taken  at  a  fixed  point  (x,y)  in  physical  space,  the  time 
derivative  on  the  right  is  at  a  fixed  point  (  ,'i)  in  the  transformed  space, 
i.e.,  at  a  given  grid  point  rather  than  at  a  given  physical  location.  Equa¬ 
tions  16-18  would  each  contain  two  additional  terms  reflecting  the  influence 
o!  the  movement  of  the  grid  points.  Thus  the  implementation  of  an  adaptive 
grid  in  the  hydrodynamic  solution  would  merely  involve  the  inclusion  of  the 
additional  terms  and  the  computation  of  ux/at  and  dy/dt  at  each  grid  point 
in  the  transformed  rectangular  plane.  These  derivatives  can  be  viewed  as  grid 
speeds  and  are  computed  from  the  new  locations  of  the  grid  points  obtained 
from  the  adaptive  grid  generator. 

31.  By  implementing  an  adaptive  grid  through  the  transformation  of  the 
time  derivative,  the  grid  and  the  flow  solution  can  be  implicitly  coupled. 

This  would  be  accomplished  by  selecting  a  grid  generation  technique  that  em¬ 
ploys  a  set  of  differential  equations  for  the  (x,y)  locations  of  the  grid 
points  which  are  dependent  upon  solution  variables,  e.g.,  the  water-surface 
elevation.  These  differential  equations  could  then  be  solved  simultaneously 
with  the  governing  flow  equations  to  yield  a  solution  in  which  the  grid  and 
the  flow  are  computed  at  the  same  time  level.  Obviously,  an  explicit  imple¬ 
mentation  could  also  be  accomplished  by  using  information  about  the  flow 
soLution  from  the  old  time  to  compute  a  new  grid.  This  approach  would  have 
to  be  taken  if  the  grid  is  not  recomputed  at  each  time-step.  As  Dwyer,  Smooke, 
and  Kee  (1982)  note,  as  the  number  of  dependent  variables  increases  and/or  the 
problem  becomes  more  nonlinear,  the  implicit  implementation  of  adaptive  grids 
becomes  less  practical  than  an  explicit  implementation. 

Fi n i le  E 1 ement  Model ing 


A  finite  element  a  1  jgorithm 

34.  Consider,  for  example,  the  depth-averaged  2-D  equations  for  free- 
surface  hydrodynamic  f 1 ows  presented  in  Appendix  A.  A  finite  element 
(FE)  algorithm  for  these  equations  formal Iv  slates  the  requirement  to  estab¬ 
lish  an  approximation  </'  (■)  to  the  toIiiI  ion  si  i  q(x,t) 


'!(*),  u  ( •  )  }  , 


to  be  constructed  on  a  discretization  Q  of  the  solution  domain 

2 

R  x  t  e{x.:  1  <_  i  <_  2,  t  ^  t  }  .  This  approximation  is  constructed  as 

0 

the  union  of  elemental  contributions  q  (x^,t)  ,  which  in  turn  are  defined 
as  expansions  in  the  k^1  degree  polynomial  cardinal  basis  set  {N^(*))  and 
time-dependent  nodal  expansion  coefficients  { Q ( t ) } ^  ,  cf.  Baker  (1983,  Ch  5). 
Mathematically,  this  statement  of  approximation  is 

qh  y  qe(xift)  (20) 

where 

q6  -  fNk(x.)}1{Q(t)}e  (21) 

35.  Since  Equations  20-21  define  the  approximation,  substitution  into 
the  flow  equations  defines  the  associated  error.  The  second  step  of  an  FE 
algorithm  is  therefore  a  formal  statement  of  constraint  on  this  error.  The 
form  of  this  constraint  which  enjoys  an  optimal  error  estimate  for  a  linear 
elliptic  problem  statement  is  the  so-called  Galerkin  statement,  which  requires 

the  distribution  of  the  error  to  be  orthogonal  to  the  function  space  {N,  (•)) 

h  ^ 

used  to  construct  q  .  For  the  more  general  problem  statement  wherein  non¬ 
smooth  solutions  to  hyperbolic  equations  are  sought,  the  Galerkin  statement  is 
typically  augmented  with  an  additional  dissipative  penalty  constraint,  cf. 
Baker  (1983,  Ch  4).  The  mathematical  statement  of  error  constraint  is  thus  of 
the  form 

/  (N  (•)}  L(qh)  +  B*  /v{ N  (-)}LC  (qh)  =  {0}  (22) 

K  K  R 

where  ,■  is  a  parameter  set  that  can  be  optimized. 

36.  Upon  evaluation  of  the  integrals  in  Equation  22,  there  results  an 
ordinary  differential  equation  system  written  on  the  time  evolution  of  the 
expansion  coefficient  set  in  the  form 

M  jt  +  (lU>1[c]  +  1 G >  +  {C})  :  (01  (23) 
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|_MJ  =  mass  matnx 
{ U  r  =  convection  velocity  field 

[cj  =  associated  (convective  derivative)  influence  on  the  dependent 


variable  set  { Q } 

{ Cl }  =  remaining  source  terms 


For  the  t ime-accurate  solution.  Equation  21  is  employed  to  evaluate  the  Taylor 


series 


{  F> 


i,!,j+i 


{ Q 1  .  -  \t  ,•  f 0 } 

J  j+. 


(01 


(24) 


Equation  24  is  a  nonlinear  algebraic  equation  system  for  0  >  0  ,  the  solution 
statement  for  which  is  cast  using  a  Newton  iteration  algorithm  in  the  form 


f'i  <^:i 


IF}? 
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(25) 


where  p  is  tlie  iteration  index.  Ihe  solution  field  is  defined  as 


1QSP  +  fSQ(P+1 


fQfj  +  l  ,x,.j  +  l  '  '  '";'j  +  l 


(26) 


md  the  Newton  Jacob ian  is  constructed  a- 
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37.  All  computational  models  of  unsteady  physical  conservation  laws 

eventually  produce  the  algebraic  equation  system  (Equation  24)  the  numerical 
solution  of  which  (Equations  25-27)  constitutes  the  heart  of  a  code.  Defin¬ 
ing  an  explicit  time  integration  procedure,  ••  0  in  Equation  24,  this  equa¬ 

tion  system  can  he  linearized,  and  the  Jacobian  (Equation  27)  replaced  by  the 
identity  matrix,  making  the  linear  algebra  solution  (Equation  25)  a  trivial 
operation.  The  "pain"  of  this  choice  is  the  restriction  to  use  of  very  small 
integration  time-steps  At  (to  bound  numerical  evolution  error).  Hence,  many 
Lime-steps  (passes  through  the  code)  are  required  to  predict  long-time  flow 
phenomena . 

38.  Conversely,  definition  of  an  implicit  integration  procedure  (A  >  0) 
significantly  relaxes  the  t  ime-st  ep/s  i  /.c  constraint  at  the  expense  of  the 


(exact)  Jacobian  (Equation  27)  becoming  a  full  nonlinear  matrix.  Hence,  in 
comparison,  forming  and  solving  Equations  24-27  is  a  much  more  computationally 
intensive  operation,  although  it  need  be  done  fewer  times  since  At  can  be 
much  larger.  Between  an  explicit  procedure,  and  use  of  the  exact  Jacobian 
with  an  implicit  procedure,  lies  the  entire  family  of  approximate  Jacobian 
methods,  e.g.,  SOR,  SLUR,  Picard,  ADI,  Checkerboard,  Zebra,  approximate  fac¬ 
torization  (AF) ,  tensor  product  (TP),  etc.  By  choosing  simple  approximations 
(block  tridiagonal),  the  degraded  (Newton)  convergence  rate  is  compensated  by 
code  operations  that  are  significantly  less  computer-resource  intensive,  and 
hence  is  (usually)  more  efficient  than  either  explicit  or  exact-implicit 
methods.  Finite  difference  code  developers  are  "expert"  at  implementing  ap¬ 
proximate  linear  algebra  solution  techniques. 

19.  Upon  substitution  of  Equation  23,  Equation  24  takes  the  specific 

form 
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The  extension  of  the  finite  element  formulational  statement.  Equations  20-22, 
to  an  adaptive  mesh  framework  requires  only  a  modest  extension  in  the  defini¬ 
tion  of  the  basis  1  (*))  upon  which  the  approximation  is  cast.  The  argu¬ 


ment  (•)  of  the  basis  involves  polynomials  in  the 


x . 


coordinate  svstem  span- 


? 


n ing  E~ 


multiplied  by  combinations  of  the  coordinates  X  of  the  nodes  of 
For  example,  in  one  dimension,  the  simplest  element 
domain  R‘  is  the  line  segment  with  a  node  at  the  left  and  right  ends.  I)e- 


the  elemental  domain. 
,1 
e 


noting  the  corresponding  values  of  X  as  X, 
the  length  of  the 
tlie  specific  form 


and  X  ,  and  realizing  that 
i,  R  j 

the  length  ot  the  element  is  A  X„  -  X,  ,  the  k=l  basis  on  R  takes 

e  K  E  e 


40.  J  n  Uk'  ox  t  oils  i  on  to  an  adaptive  grid,  the  coordinates  of  the  nodes 
ot  the  nu-sh  are  permitted  (required)  to  move  (in  time),  hence  X  -  X(t)  . 
Therefore  Equation  21  is  generalized  to  the  form 


‘ll  (*,t)  1  Nk(x,X(t))  ;  1  lQ(t  ) 


(30) 


Substitution  of  Equations  20  and  30  into  Equation  22  and  using  the  chain  rule 
for  the  time  derivative  yields 
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where 


svstevn.  for  e 


is  a 
xamp  1  e 


matrix  with  elements 
for  Equation  20 


a 


function  of  the 


x  coordinate 


(32) 


where  the  elements  of  X-  express  the  (Eagrangian)  velocity  of  the  nodal 
coordinates  of  the  adaptive  grid.  1'pon  substitution  of  Equation  31  into  the 
basic  finite  element  constraint  statement.  Equation  22, 


and  proceeding  through 


Che  formulational  details,  for  implementation  of  an  adaptive  grid  Equation  28 
becomes  modified  to  the  form 
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1'hus  the  motion  of  the  adaptive  grid  is  accounted  for  by  tile  grid  velocity 
action  on  the  solution  variable  lOl  through  the  grid  convection  matrix 

M  • 

41.  The  step  to  use  of  a  boundary-fitted  grid  (BFG)  transformation,  as 
described  in  detail  in  Baker  (1983,  Cli  3,8),  involves  transformation  of  the 
divergence  operator  into  scalar  components  parallel  to  the  coordinate  lines 
of  the  BFC;  system.  The  associated  scalar  product  with  the  convection  velocity 


vector  u.  produces  scalar  components  that  are  also  parallel  to  the  BFG  co¬ 


ordinate  lines.  Denoting  this  system  as  ■  ,  and  proceeding  through  the  de¬ 

tails  produces  Equation  33  in  the  modified  form 
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In  Equation  34,  1 1  jo  discrete  index  K  signifies  that  all  derivatives  and  con¬ 
vection  velocity  components  are  computed  in  scalar  components  parallel  to  n 
(rather  than  >;  ^ )  . 

Gomputat ional  aspects 

42.  Implementation  of  the  finite  element  algorithm  statement,  Equa- 
t  ions  20-22  and  29-27,  using  the  adaptive  gr id  generalization  and  the  BFG 
transformation  provides  for  several  computational  aspects  of  importance.  The 
inclusion  of  the  grid  convert  ion  contrihut  ion,  in  Equations  33-34  ,  per¬ 

mits  ile  f i n i t  ion  and  use  of  adapt  ive  methodologies  to  cluster  refined  grids 


about  iiu  v  i  iii;  t  routs.  A  key  implementation  aspect  relates  to  code  efficiency. 

I  mite  i  lenient  I  laid  mechanics  codes  have  gained  a  reputation  for  being  less 
ell  i  c i on  L  than  their  finite  difference  counterparts,  As  noted  earlier,  this 

.  I  I  1  'jtieiice  ol  1  lu  llnite  ,•  |<  neilt  theory  itself,  hut  instead  reflects 
directly  on  implementation  choices  for  the  linear  algebra  statement,  Equa¬ 
tions  2b-27  fur  an  implicit  construction,  ■  '•  0  in  Equation  24,  as  discussed 
Based  on  the  finite  element  structural  mechanics  historical  development  of 
direct,  out-of-core  solvers,  many  (most)  finite  element  code  implementations 
nave  i\  I  ied  on  direct  extension  of  these  solver  software  packages  to  the  fluid 
mechanics  problem  classes.  In  distinct  ion,  a  direct  solution  (of  even  a 
linear  Poisson  equation)  is  never  attempted  in  an  FD  code,  but  instead  a  ma¬ 
tt  i  iteration  procedure  is  defined  that  uses  easily  formed  (usuallv  block- 
tridiagonal)  approximations  to  the  Newton  algorithm  Jacobian  (Equation  27). 

I  lie  computer  storage  requirement  for  the  approximate  Jacobian  is  negligible, 
in  comparison,  as  is  the  CPU  needed  to  execute  an  LU  decomposition  and  hack 
sub.-;  t  i  tu  t  i  on  .  Hie  computational  penalty  is  a  much  reduced  convergence  rate 
for  Equation  2b,  but  each  iteration  proceeds  so  rapidly  that  the  many  grid 
sweeps  required  for  the  El)  code  are  typical]'.'  completed  with  minimal  expendi¬ 
ture  of  computer  resources. 

-Vi.  The  generalized  (BEG)  coordinates  form  of  the  finite  element  al¬ 
gorithm,  Equation  14,  provides  the  theoretical  framework  for  implementation  of 
efficient  solving  techniques  for  comparison  to  El)  procedures.  In  particular, 
as  developed  in  detail  by  Baker  ( 1 1)8  ) ,  CL  b)  for  the  BEG  form  of  the  algorithm 
statement,  using  the  tensor  product  cardinal  basis  fN^(*)t  for  Equation  21, 
hence  i  » ,  facilitates  reexpression  of  the  Newton  algorithm  Jacobian  £0  > 
t  ion  27,  in  the  outer  (tensor)  matrix  product  form 

[|]  [>,]  g  [|J  on 


Then,  for  example,  for  file  lie.  ai  basis  ( k—  l )  and  the  governing  vertically 
averaged  hvdrodyn.imic  equal  ion-.,  each  tensor  product  factor  [']  in  Equa- 
t  ion  3b  is  a  three-block  tridiagonal  matrix,  the  storage  requirements  for 
which  are  ncgl  irihle  in  ■  ■  to  .i  ■  i  ■  •  lo  that  of  [~f]  ,  as  is  the  CPU.  execution 
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which  is  replaced  in  the  code  by  the  approximation 
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Equation  37a  is  solved  by  sweeping  the  grid  lines  parallel  to  the  (BFG) 

coordinate  to  obtain  the  intermediate  vector  {P}  ,  the  entries  in  which  are 
then  realigned  parallel  with  the  n,;  coordinate  for  subsequent  solving  of 
Equation  37b.  This. second  solution  sweep  yields  the  tensor  product  approxima¬ 
tion  to  { 5Q}  j+j  »  which  is  then  used  in  Equation  26  to  compute  {Q}^j  .  The 
Newton  algorithm  quadratic  convergence  rate  is  compromised,  of  course,  by  using 
the  approximation  defined  by  Equation  37,  but  it  is  (usually)  more  than  com¬ 
pensated  by  the  much  reduced  computer  resource  requirements. 

45.  The  definition  and  use  of  the  tensor  matrix  product  approximate 
Jacobian  exert  another  influence  of  great  importance.  With  elimination  of  the 
Newton  Jacobian  as  the  dominating  computer  (hardware)  impediment,  a  finite 
element  code  can  then  afford  to  use  a  much  finer  discretization  of  the  hy¬ 
draulics  problem  domain  along  with  an  adaptive  grid.  The  net  result  would  be 
a  significant  improvement  in  solution  accuracy  with  enhanced  definition  of 
sharp  field  gradients.  This  in  turn  would  most  likely  decrease  sharply  the 
need  for  excessive  eddy  (or  artificial)  diffusivity,  yielding  an  overall  sub¬ 
stantial  accuracy  improvement. 
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PART  HI:  GENERATION  OF  ADAPTIVE  GRIDS 


46.  The  generation  of  time-dependent  grids  to  better  resolve  developing 
; low  gradients  can  be  accomplished  while  either  maintaining  a  fixed  number  of 
grid  points  or  adding  grid  points  in  locations  where  some  algorithm  that  senses 
local  error  indicates  they  are  needed.  Most  of  the  research  over  the  past  few 
years  on  the  use  of  adaptive  grids  has  adopted  the  fixed  number  of  points  ap¬ 
proach.  For  completeness,  a  brief  discussion  of  the  variable  number  of  points 
approach  is  also  given. 

47.  In  order  to  construct  a  suitable  grid  to  be  used  for  solving  a 
particular  problem,  some  knowledge  of  the  solution  is  required  initially  if 
the  grid  is  to  remain  fixed  in  time.  For  example,  a  high  grid-point  density 
is  desired  in  high  gradient  regions  of  the  solution  variables.  However,  in 
general,  too  little  is  known  about  the  solution  to  be  computed;  and  the  loca¬ 
tions  of  such  high  gradient  regions  are  unknown .  Therefore  the  establishment 
of  a  suitable  grid  a  priori  is  difficult. 

Fixed^  Number  of  Gri_d  Points 

48.  Accepting  that  the  best  grid  for  a  given  problem  depends  on  the 
solution,  a  variety  of  papers  concerned  with  moving  the  grid  points  as  the 
solution  changes  have  been  published  by  researchers  in  the  aerodynamics  field. 

A  basic  review  paper  by  Thompson  (1984)  presents  an  extensive  discussion  of 
these  publications.  In  addition,  basic  concepts  of  grid  generation  techniques 
are  presented  by  Thompson  (1983c).  'The  purpose  of  the  present  discussion  is 
not  to  duplicate  Thompson's  review  but  rather  to  focus  on  the  two  basic  ap¬ 
proaches  for  generating  adaptive  grids. 

Va r  i a  t  i  ona  1  a pproac. It 

49.  A  general  feature  of  all  adaptive  grid  generation  techniques  is  the 
minimization  of  numerical  error  bv  the  equ idistribut ion  of  some  solution 
quantity  related  to  the  error  over  the  grid.  From  a  continuous  viewpoint, 
since  some  quantity  is  to  be  minimized  by  the  movement  of  the  grid,  a  varia¬ 
tional  approach  appears  logical.  As  will  be  shown  later,  the  grid  generation 
approach  developed  by  Thompson  i-t  al.  (1974),  is  actually  a  special  case  of 
grid  generation  based  upon  a  variational  principle. 


30.  NmiRTii.il  iTior.  Hoi  or  l.-  di  .cuss  i  ng  some  specifics  of  grid  genera¬ 
tion,  a  discussion  of  solution  errors  due  to  the  grid  is  required.  As  demon¬ 
strated  h v  Mast  in  (1M82),  and  Thompson  and  Mastin  (1983),  a  curvilinear 
coordinate  system  can  have  a  significant  effect,  favorable  or  adverse,  on  the 
error  in  the  numerical  solution.  This  is  because  the  truncation  error  is  de¬ 
pendent  not  only  on  the  higher  order  derivatives  of  the  solution  and  the  local 
grill  spacing  but  also  on  the  rate  of  change  of  the  grid  spacing  and  on  the 
skewness  of  the  grid,  i.e.,  the  departure  of  the  grid  lines  from  orthogonality. 

31.  In  Mast  in's  analysis  it  is  shown  that  using  standard  central  dif¬ 
ferences  to  approximate  first  derivatives  in  the  computational  plane  results 
in  additional  truncation  errors,  such  as  that  given  below  for  f.  ,  due  to  the 
nonuniformity  and  nonorthogonality  of  the  grid. 

T  =  OOO  -  1/2  x _ x _  .  f  -  1/2  (x.y..  +  y.x.  .)  f 

xx  ■  \  v.  xy 

-  1/2  y.v. _f  +  H.  0.  T.  (38) 

yy 

where 

h  =  the  grid  spacing 
f  =  an  arbitrary  variable 

Considering  a  rectangular  grid,  if  the  grid  spacing  in  the  x  direction  is 

uniform  then  x..  =  0  and  furthermore  since  the  grid  is  orthogonal  then 

y.  =0  ,  which  results  in  the  vanishing  of  the  additional  terms  above  and  thus 

2 

tile  truncation  error  being  0(h  )  .  However,  in  general  these  additional  terms 
do  not  vanish  but  can  be  shown  to  he  proportional  to  the  rate  of  change  of 
the  grid  spacing  and  to  the  departure  from  orthogonality.  Note  that  since  the 
additional  terms  involve  second  derivatives  of  the  solution  variable,  they 
constitute,  in  effect,  a  negative  numerical  diffusion.  It  is  observed  that  un¬ 
less  care  is  exercised  in  generating  the  coordinate  system,  difference  methods 
for  first-order  partial  differential  equations  will  be  only  first  order  and 

those  for  second-order  equations  will  he  inconsistent.  However,  the  trunca- 

2 

tion  error  for  first-  and  second-order  derivatives  can  be  increased  to  0 (h"“ ) 

and  0(h)  ,  respectively,  if  the  second-order  derivatives  of  x  and  v  are 
> 

0(h)  .  Ibis,  of  course,  I imits  the  rate  of  change  in  grid  spacing  and  the 
curvature  of  the  coordinate  lines. 

32.  Mastin  further  shows  that  the  degree  to  which  nonorthogonality 


increases  the  truncation  error  is  determined  by  the  value  of  the  Jacobian  of 
the  mapping.  The  factor  that  causes  an  increase  in  truncation  error  for  first 
derivatives  is  derived  to  be  1/sin  ("  -  -f)  where  i)  and  0  denote  the 
angles  of  inclination  of  the  o  =  constant  and  '  =  constant  coordinate  lines. 
Ibis  is  exactly  the  factor  by  which  the  Jacobian  is  decreased  by  the  shearing 
of  the  coordinate  lines.  A  reasonable  degree  of  nonorthogonality  has  a  neg¬ 
ligible  effect;  however,  the  rate  of  increase  in  truncation  error  increases  as 
the  grid  becomes  less  orthogonal.  It  can  be  shown  that  when  the  angle  between 
grid  lines  is  greater  than  45  deg  the  error  from  nonorthogonality  is  less  than 
that  from  a  nonuniform  grid  spacing.  In  addition,  computations  have  demon¬ 
strated  that  doubling  the  grid  spacing  from  one  point  to  the  next  results  in  a 
17  percent  error  in  flow  solutions.  Thus  the  rate  of  change  of  the  grid  spac¬ 
ing  and  the  departure  from  orthogonality  must  be  controlled,  or  else  signifi¬ 
cant  error  can  be  introduced. 

53.  Constraints ■  The  above  discussion  on  grid-induced  errors  has  been 
given  as  a  lead-in  to  the  generation  of  adaptive  grids  using  a  variational 
approach.  On  the  one  hand,  there  is  a  need  to  force  grid  points  to  concentrate 
near  large  solution  gradients;  on  the  other  hand,  as  illustrated  in  the  previ¬ 
ous  discussion,  there  is  a  need  to  generate  grids  that  are  relatively  smooth 
and  don't  deviate  too  much  from  being  orthogonal.  The  calculus  of  variations 
is  well-suited  to  handle  such  problems  since  integrals  over  the  grid  can  be 
written  which  measure  the  three  desired  features  discussed  above,  namely, 

(a)  concentration  of  grid  points  near  large  solution  gradients,  (b)  a  smooth 
distribution  of  grid  points,  and  (c)  a  relatively  orthogonal  grid.  The  varia¬ 
tional  formulation  of  adaptive  grids  is  discussed  in  detail  by  Brackbill 
(1982).  A  brief  discussion  follows  here. 

54.  The  volume  of  computational  cells  in  three  dimensions  (the  area  in 
two  dimensions)  is  given  by  the  Jacobian,  J  ,  of  the  mapping  (Equation  6). 
Therefore,  if  the  integral 

]  =  f w(x,v)  J  dV  ,  (39) 

w  J 

I) 

which  is  a  measure  of  the  weighted  variation  of  the  computational  cell  size 
over  the  grid,  is  minimized  for  some  weight  function  w(x,v)  ,  which  is  some 
measure  of  solution  error,  a  concentration  of  grid  points  can  be  obtained.  In 
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other  words,  where  the  weight  function  w(x,y)  becomes  large  the  cell  size 
becomes  small,  and  where  w(x,y)  becomes  small  the  size  of  the  computational 
cell  becomes  large. 

55.  Likewise,  the  smoothness  of  the  grid  is  measured  by  the  integral 


(40) 


with  the  orthogonality  measured  by 

,  ■  /<• 

I) 


2  3 

Vri)  J  dV 


(41) 


where  the  J  is  added  to  cause  orthogonal ity  to  be  emphasized  more  in  larger 
cells.  Note  that  \j\  •  Vri  =  0  for  an  orthogonal  grid.  Therefore  an  adaptive 
grid  generator  can  be  developed  by  minimizing  a  sum  of  the  three  integrals 
given  in  Equations  39-41,  i.e.. 


I  =  A  I  +  A  I  +  A  I 
s  s  o  o  w  w 


56.  From  the  above  expressions  for  1,1,  and  I  we  have 

sow 


2  5  3 
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where 

N  =  a  characteristic  number  of  points 
L  =  a  characteristic  length 
W  =  average  weight  function  over  the  field 
Thus,  for  Equation  42  to  be  dimensionally  correct 


where  A'  and  A'  are  positive  constants.  In  two  dimensions,  the  factors 
7  °  W5  4 

(N/l.)  and  (N/L)  both  become  (N/L)  since  the  Jacobian  is  then  propor- 
2  2  3  3 

tional  to  l.^/N  rather  than  to  L  /N  ,  The  characteristic  length  and  number 
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ot  [mint.-?  m  L  >.;h  t  logically  be  taken  as  the  cube  roots  of  the  volume  and  the 
total  number  of  points  in  the  field,  respectively,  in  three  dimensions,  the 
square  root  being  used  in  two  dimensions. 

.  obviously,  bv  selecting  appropriate  values  of  X  ,  X'  ,  and  X'  , 

s  w  o 

grids  that  emphasize  smoothness,  concentration  of  grid  points,  or  orthogonal¬ 
ity  can  be  generated.  The  variational  formulation  thus  provides  a  competitive 
enhancement  of  these  three  grid  characteristics.  A  note  of  caution  concern¬ 
ing  orthogonality  is  perhaps  needed.  Purely  orthogonal  grids  cannot  be  gen¬ 
erated  '.'liiMi  prescribing  the  location  of  boundary  points,  even  if  the  condition 
•  ■  =  0  is  enforced.  Since  derivatives  of  the  coordinates  have  to  be 

spec i tied  to  satisfy  orthogonality  at  the  boundaries,  specifying  the  location 
of  the  points  overspecifies  the  problem  with  second-order  partial  differential 
equations  as  the  grid  generation  system.  Therefore,  in  order  to  generate  a 
strictly  orthogonal  grid,  the  boundary  points  must  be  allowed  to  move  on  the 
boundary.  Obviously,  in  many  cases  this  constitutes  a  major  restriction  on 
the  generation  of  a  useful  grid.  Therefore  large  values  of  X^  compared  with 

values  of  X  and  > ’  will  only  enhance  the  orthogonality  of  the  grid.  In 
s  w 

general,  the  grid  is  still  nonorthogonal  and  all  terms  in  the  governing  trans¬ 
formed  partial  differential  equations  must  be  retained. 

58.  For  a  2-1)  adaptive  grid  generator,  partial  differential  equations 
for  the  solution  of  the  physical  coordinates  (x,y)  are  desired.  These  are 
obtained  through  the  calculus  of  variations  applied  to  the  minimization  of  the 
sum  of  integrals  in  Equation  42.  In  general,  if  we  wish  to  find  functions 
:.(x,y)  that  are  differentiable  on  (x,y)  and  satisfy  fixed  constraints  on 
the  boundary  of  the  domain  that  minimize  some  integral  functional 


I) 


the  calculus  of  variations  gives  that  !.(x,y)  are  the  solutions  of  the  Euler 
Lagrange  equations 


(44) 


59.  Two-dimensional  generator.  liraokbill  (1982)  has  derived  the  2-D 


Euler-Lagrange  equations  for  the  minimization  of  Equation  42  as 


b.x_  +  b.,x.  +  b..x  4-  a .  v .  .  +  a..y  +  a^y 
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+  anx 
3  nn 


+  clV; 


+  c  y.  +  c_y 

2  .\n  3'nn 


2  Aw  3w 
J  2  3y 


For  completeness,  the  coef f icients  are  reproduced  from  the  cited  paper  and 
presented  below: 


a  =  \  a  +  A  wa  +  A  a 
i  s  s  .  w  w .  o  o  . 

ill 


b .  ■  \  b  +  A  vb  +  A  b  ,  and 
r  s  s  .  w  w .  o  o . 

l  l  i 


c .  =  A  c  +  wc  +  A  c 
l  s  s  .  w  w .  o  o  . 

l  i  l 
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=  -At 

5  bsl 
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Csl 

=  C  t 

as2 

=  2  a;-; 

;  bs2 

=  -2BH  ; 

Cs2 

=  -2C 

as3 

=  -Ay 

;  bs3 
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°s3 
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where 


A  =  x .  y .  +  x  y 
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and 


(46) 
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a  =  x . y  +  x  y.  ,  b  „  =  -2y  v  ,  c  0  =  -2xrx 

\.l )  n  re  7  r  n  u/  r 


n  n 


■  .  n 


aw3  =  "X-y- 


b  _  =  v. 
w3 


c  _  =  x . 
w3  >. 


a  -  =  x  y 
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b  =  x 
ol  M 


col  =  yn 


a  „  =  x_y  +  x  y.  ,  b  =  2(2xrx  +  y  v  )  ,  c  „  =  2(x^x  4-  2y^y  ) 

o2  /n  n  .  o2  r  n  • .  n  °2  r.  n 


o3  "  ‘ 


x.y. 


b  -  =  x . 
o3 


co3  =  yf; 


60.  Equation  45  constitutes  the  adaptive  grid  generator  in  2-D.  Finite 
differences  could  be  used  to  solve  the  equations  on  the  transformed  rectangular 
plane  for  values  of  (x,y)  at  each  mesh  point.  These  of  course  are  then  used 

to  compute  the  metrics,  e.g.,  x.  ,  y^  ,  etc.,  in  the  transformed  differential 
equations  (Equations  16-18).  The  solution  of  Equation  45  could  be  obtained  at 
each  time-step  in  an  implicit  fashion  with  the  flow  equations  or,  if  desired, 
the  grid  generation  could  be  lagged  one  or  more  time-steps  behind  the  flow 
computations  to  yield  an  explicit  grid  generator.  In  addition,  as  previously 
discussed,  interpolation  of  solution  variables  from  the  old  grid  to  the  new 
grid  could  be  used  to  implement  the  adaptive  grid  or  the  implementation  could 
be  accomplished  by  transforming  the  time  derivative  in  the  governing  equations. 

61.  An  interesting  development  occurs  if  the  coefficients  in  Equation  42 


are  set  to  be 


X  =  1 
s 

X’  =  X'  =  0 
w  o 


(47) 


to  generate  as  smooth  a  coordinate  system  as  possible.  If  the  algebra  is  per¬ 
formed,  it  is  seen  that  the  grid  generating  system  given  by  Equation  45  becomes 

ax. .  -  2 Ex  +  yx  =0 

•  n  nn 

(4 
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where 


2.2 

i  =  x  +  v 

n  ■  n 

.  j  =  x . x  +  v  _ v  (49) 

'I  -  */ n 

2  1 

i  =  xT  +  yZ 

In  other  words,  the  generating  system  becomes  the  elliptic  system  given  by 

Equations  5  and  b  where  the  P  and  <)  control  functions  are  set  to  zero, 

i.e.,  an  elliptic  system  based  upon  Laplace's  equation  rather  than  Poisson's 

equation  is  used.  Similarly,  if  i  0  the  weighting  function  w(x,y)  , 

w 

which  as  noted  is  related  to  some  measure  of  the  solution  error,  can  be  re¬ 
lated  to  the  control  functions  in  an  elliptic  system  based  on  Poisson's 
equat  ion . 

bd.  One-d  iinens  imia  1  generator.  An  adaptive  grid  generator  based  upon 
the  variational  approach  can  also  be  derived  for  1-D  problems,  e.g.,  a  surge 
moving  down  a  river.  If  the  coefficients  in  Equation  45  are  evaluated  for  the 
case  of  x  being  the  only  dependent  variable,  the  grid  generation  equation 
becomes 


t 

s 


+ 


_w  dw 
2  dx 


(50) 


where  the  Jacobian,  J  ,  reduces  to  x .  .  Neglecting  the  smoothness  contribu¬ 
tion,  this  equation  can  be  reduced  to 


or 


d  (x_  ) 
d’ 


d_x  1  dw 
Xr  d r  2w  dx 


(51) 


d  (xr  ) 


1  dw 

2  w 


(52) 


Integrating  Equation  52  then  yields 


1 

=  -  •!!  w  +  constant 
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■  n  x  . 


(53) 


constant 


(54) 


■’n 


which  implies 


1/2 


x  _w 


constant 


(55) 


This  is  similar  to  the  generating  system  below  for  the  equidistribution  of  w 
discussed  by  Thompson  (1983b) 


x.w(x)  =  constant  (56) 

Therefore  minimization  of  error  by  equidistribution  is  equivalent  to  solving 
a  Eu  1 er-Lagrange  equation  arising  from  a  variational  approach. 

63.  Weight  functions.  As  noted,  the  effect  of  the  weight  function,  w  , 
is  to  reduce  the  grid  spacing  where  w  is  large,  and  therefore  the  weight 
function  should  be  set  as  some  measure  of  the  solution  error,  or  as  some  mea¬ 
sure  of  the  solution  variation.  The  discussion  to  follow  has  been  taken  from 
Thompson,  Warsi,  and  Mastin  (1984).  The  simplest  choice  is  just  the  solution 
gradient,  e.g.,  in  one  dimension 


w  =  u  (57) 

x 


where  u  is  some  solution  variable.  In  this  case.  Equation  56  becomes 


x.u  =  constan' 
x 


which  then  reduces  to 


u .  =  constant 


With  the  solution  gradient  as  the  weight  function,  the  point  distribution  thus 
adjusts  so  that  the  same  change  in  the  solution  occurs  over  each  grid  interval, 
as  illustrated  below: 


n 


I'li  is  choice  tor  the  weight  function  has  the  disadvantage  of  making  the  spacing 
infinitely  large  where  the  solution  is  flat,  however. 

64.  A  closely  related  choice,  also  based  on  the  solution  gradient,  is 
the  form 

(58) 

Now  an  increment  of  arc  length,  ds  ,  on  the  solution  curve,  u(x)  ,  is  given 

by 

9  9  2  2  2 

ds“  =  dx  +  du  =  (1  +  u~)  dx 

so  that  this  form  of  the  weight  function  may  be  written 

w  -  s 

x 

and  then  Equation  56  becomes 

x  s  =  constant 
x 

which  reduces  to 

s  -  cons t  an t 


Thus,  with  the  weight  function  defined  hv  Equation  58,  the  grid  point  di  -tri- 


occurs  over  each  grid  interval, 
lowing  point  distribution: 


For  the  curve  shown  above  this  gives  the  fol- 


x 


Unlike  the  previous  choice,  this  weight  function  gives  uniform  spacing  when 
the  solution  is  flat.  The  concentration  of  points  in  the  high-gradient  region, 
however,  is  not  as  great. 

65.  This  concentrat ion  can  be  increased,  while  still  maintaining  uniform 
spacing  where  the  solution  is  flat,  hv  altering  the  weight  function  to 


w 


•>  a 

<~u‘" 

x 


(59) 


where 

against 


is  a  parameter  to  he  specified.  Mow  considering  u  to  be  plotted 
x/ t  ,  we  have  for  an  increment  of  arc  length  on  this  solution  curve, 


so  that  this  weight  function  is  equivalent  to 


w  -  s 


and  Kquat ion  56  becomes 


constant 


whii'li  rcdiu'i'K  to 


s  _  =  constant 

Thus  wo  have  equal  increments  of  arc  length  on  the  solution  curve  with  u 

plotted  against  x/ <  in  this  case.  Now  division  of  the  abscissa  by  t  for 

a  flat  curve  would  simply  reduce  the  spacing  by  the  same  factor.  However, 

since  the  slope  steepens  as  the  curve  is  compressed  to  the  left  by  this  change 

of  scale,  the  effect  on  the  spacing  where  the  curve  is  not  flat  will  be  a 

greater  reduction  in  spacing.  In  fact,  since  the  1  in  the  weight  function 

2  2 

given  hv  liquation  59  tends  to  produce  equal  spacing,  while  the  x  u  tends 
to  produce  concentration  in  the  high-gradient  regions,  witli  infinite  spacing 
in  I  lat  regions,  this  weight  function  involves  a  weighted  average  between  the 
tendency  toward  equal  spacing  and  that  toward  concentration  entirely  in  the 
high-gradient  regions.  The  larger  the  value  of  t  ,  the  stronger  will  be  the 
concent  rat  ion  in  the  high-gradient  regions  and  the  wider  the  spacing  in  the 
flat  regions. 

hh .  A  disadvantage  of  all  the  above  forms  of  the  weight  function  is 
that  in  regions  near  solution  extrema,  i.e.,  where  =  0  locally,  spacing 

is  similar  to  ! hat  in  flat  regions,  as  is  illustrated  below  for  the  form  given 
by  F.qua  t ion  57: 


x 


Although  Liu*  d  is  t  r  ihu  t  ions  produced  by  the  solution  arc  length  forms,  liqua¬ 
tions  5>H  and  would  have  closer  spacings  near  the  extrema,  the  effect  is 

still  r  1  i « ■  same,  i.e.,  to  concentrate  points  on  1 v  near  gradients,  not  extrema. 


67.  Concentration  near  solution  extrema  can  be  achieved  by  incorporating 
some  effect  of  the  second  derivative,  ,  into  the  weight  function.  A  logi¬ 

cal  approach  is  to  include  this  effect  through  consideration  of  the  curvature 
of  the  solution  curve: 


K  =  -- 


(' +  u7 


If  the  weight  function  is  taken  as 

>  i  I 

w  =  1  +  a“|K  |  (60) 

points  will  be  concentrated  in  regions  of  high  curvature  of  the  solution  curve, 
e.g.,  near  extrema,  with  a  tendency  toward  equal  spacing  in  regions  of  zero 
curvature,  i.e.,  where  the  solution  curve  is  straight  (not  necessarily  flat). 
This  weight  function,  however,  has  the  serious  disadvantage  of  treating  high- 
gradient  regions  with  little  curvature  essentially  the  same  as  regions  where 
the  curve  is  flat.  Thus  in  the  curve  shown  above,  nearly  all  the  points  would 
be  concentrated  near  the  maximum  in  the  curve,  with  very  wide  spacing  in  the 
high-gradient  regions  on  both  sides. 

68.  A  combination  of  the  weight  functions  given  by  Equations  59  and  60 
provides  the  desired  tendency  toward  concentration  both  in  regions  of  high 
gradient  and  near  extrema: 

w  =  (1  +  li2\K\)yjl  +  u2u^  (61) 

where  i  and  8  are  parameters  to  be  specified.  Clearly,  concentration  near 
high  gradients  is  emphasized  bv  large  values  of  n  ,  while  concentration  near 
extrema  (or  other  regions  of  large  curvature)  is  emphasized  by  large  6  . 

69.  Another  approach  to  the  inclusion  of  the  second  derivative  is 
simpLy  to  take  the  weight  function  as 


u 

+  6  u 

X 

XX 

where  i  and  8  are  nonnegative  parameters  to  hi1  specified.  Note  that  the 
replacement  of  Equal  ion  5‘J  with  the  form  given  by  Equation  62,  with  8  -  0  , 


still  leaves  a  reasonable  form  for  the  weight  function,  but  the  clear  associa¬ 
tion  with  the  geometric  properties  of  the  solution  curve  is  lost.  In  this 
case,  the  weight  function  corresponding  to  Equation  58  would,  after  substitu¬ 
tion  in  Equation  56,  lead  to  the  condition 

x,  +  u„  =  constant 

which  corresponds  to  an  equal  distribution  of  the  distance  between  points  on 
the  solution  curve  along  a  right-angle  path  formed  by  Ax  and  Au  from  one 
point  to  the  next.  While  this  distance  has  some  indirect  relation  to  arc 
length  on  the  solution  curve  (the  chord  length  being  the  hypotenuse  of  the 
right  triangle  formed  by  this  Ax  and  Au),  the  direct  association  with  arc 
length  is  preferable.  Following  the  same  reasoning,  the  use  of  solution  curve 
curvature,  rather  than  simply  the  second  derivative,  is  also  preferable. 
Therefore  the  form  given  by  Equation  61  is  probably  more  appropriate  than 
that  of  Equation  62. 

70.  Since  the  numerical  evaluation  of  higher  derivatives  can  be  subject 
to  considerable  computational  noise,  the  use  of  formal  truncation  error  ex¬ 
pressions  as  the  weight  function  is  usually  not  practical,  hence  the  emphasis 
above  on  solution  gradients  and  curvature.  Some  problems  may  arise  even  with 
solution  curvature,  i.e.,  with  second  derivatives,  in  rough  transits.  It  is 
common  in  any  case  to  limit  the  grid  point  movement  at  each  time-step  and/or 
to  smooth  the  new  point  distribution. 

71.  One  final  note  on  weighting  functions  might  be  offered.  When  more 
than  one  dependent  solution  variable  is  being  computed,  e.g.,  the  water-surface 
elevation  and  u  and  v  velocities  in  depth-averaged  flow  computations,  it 

is  not  clear  how  to  choose  the  most  appropriate  weight  function.  Some  re¬ 
searchers  have  generated  separate  adaptive  grids  for  solution  variables,  e.g., 
Dwyer,  Smooke,  and  Kee  (1982),  whereas  others  such  as  Yanenko  et  al .  (1982) 
have  attempted  to  use  weight  functions  that  are  composed  of  a  dependence  on 
several  solution  variables  simultaneously.  Additional  comments  on  weight  func¬ 
tions  for  free-surface  hydrodynamic  problems  will  be  offered  in  PART  IV. 

A t tract i on- repul s ion  ap p r o  a c h 

72.  In  this  approach  to  adaptive  grid  generation,  the  grid  points  move 
directly  under  the  influence  of  mutual  attraction  or  repulsion  between  points. 
The  basic  idea  is  to  generate  a  grid  speed  for  each  grid  point  as  a  result  of 


38 


the  collective  attraction,  or  repulsion,  of  all  other  points.  This  is  ac¬ 
complished  by  assigning  to  each  point  an  attraction  (repulsion)  that  is  pro¬ 
portional  to  the  difference  between  some  measure  of  the  local  truncation  error 
uid  a  measure  of  the  average  truncation  error  over  the  complete  grid.  Thus 
points  with  a  local  error  that  exceeds  the  average  error  attract  other  points, 
with  a  consequent  reduction  in  grid  spacing  and  thus  a  reduction  in  local  er¬ 
ror,  while  points  with  a  local  error  less  than  the  average  repel  other  points 
and  thus  increase  the  local  grid  spacing.  In  analogy  with  the  electrostatic 
force  between  two  electrical  charges,  the  attraction  exerted  by  a  point  on 
other  points  is  attenuated  by  an  inverse  power  of  the  distance  between  points 
in  the  computational  region. 

73.  One-d imensional .  For  two  points  A  and  B,  Anderson  and  Rai  (1982) 
write  the  grid  velocity  at  B  due  to  error  at  point  A  in  computational  space 


vba  =  k 


S|a  -  I  6  lav 


where 


K  =  a  proportionality  factor  chosen  so  that  the  velocity  at  any  point 
does  not  exceed  a  preset  maximum 

?|  =  a  measure  of  the  truncation  error 


With  V  =  ^  ,  the  grid  speed  in  1-D  physical  space  becomes 


If  the  number  of  fixed  points  is  N  ,  a  summation  over  all  points  yields  the 
value  of  x^  to  be  inserted  into  the  transformed  flow  equations  with  the  time 
derivative  transformed  bv  Equation  19.  A  time  integration  of  the  grid  speed 
equation  then  yields  the  physical  location  of  the  new  grid  points  for  use  in 
computing  the  metrics  of  the  transformation. 

74.  Two-dimensional.  A  dismission  of  a  2-1)  adaptive  grid  generator 
based  upon  the  concept  of  point  attraction  or  repulsion  is  also  presented  by 
Anderson  and  Rai  (l')87).  V,  in  tin  1  M  ,  a-.e,  .aid  -.poods  and  v  are 


determined  bv  prescribing  a  certain 


.act  ion  that  is  a  function 
■t  .  t  i  .>r  est  i mates  in  the 


and  n  directions.  A  summation  over  all  points  then  yields  the  final  grid 
speed  equations. 

75.  Error  forms.  As  in  the  variational  approach,  some  estimate  or 
measure  of  the  truncation  error  must  be  selected  to  cause  a  motion  of  the  grid. 
The  term  |e|  in  Equation  63  could,  of  course,  be  the  actual  truncation  error; 
however,  any  other  quantity  could  also  be  used.  For  example,  j  e|  could  be 
taken  to  be  any  of  the  forms  presented  for  the  weighting  function,  w  ,  in 
Equations  57-62.  The  same  comments  previously  made  concerning  the  use  of 
higher  derivatives  of  solution  variables  also  apply.  Therefore,  if  higher 
derivatives  are  used  in  the  computation  of  |e|  ,  large  amounts  of  smoothing 
will  probably  be  required. 

76.  Constraints .  A  disadvantage  of  using  adaptive  grid  generators  based 
upon  the  grid  speed  approach  is  that  control  of  grid  motion  is  not  easily  main¬ 
tained.  Such  control  is  one  of  the  major  advantages  of  using  the  variational 
approach  since  control  is  inherently  achieved  through  the  specification  of 

X^  ,  X^  ,  and  X^  in  Equation  42.  Although  grid  control  is  difficult,  the 
point  attraction  approach  will  always  prohibit  the  crossing  of  grid  lines  of 
the  same  family.  This  is  illustrated  by  Equation  63.  As  two  points  approach 
each  other,  the  local  error  becomes  less  than  the  average  error  and  thus  the 
force  that  drives  the  grid  becomes  negative,  i.e.,  repulsive,  resulting  in 
the  points  moving  apart.  Attempts  at  controlling  grid  oscillations  and  non¬ 
smoothness  of  the  grid  have  been  made  by  damping  the  grid  speeds  permitted  by 
limiting  the  change  allowed  in  the  Jacobian,  i.e.,  cell  size,  at  each  point. 
This  technique  is  discussed  in  detail  by  Anderson  and  Rai  (1982). 

Variable  Number  of  Grid  Points 

77.  Most  of  the  research  in  adaptive  grids  has  followed  the  approach  of 
moving  a  fixed  number  of  grid  points  in  order  to  be  aBle  to  accurately  resolve 
steep  fronts  in  flow  computations,  or  perhaps  concentration  profiles,  during 
transport  processes.  This  approach  utilizes  coordinate  transformations  so 
that  even  though  the  grid  points  are  moving  in  the  physical  plane  all  compu¬ 
tations  are  performed  in  a  fixed,  rectangular,  transformed  plane.  Another 
approach  is  to  add  or  subtract  grid  points  as  the  physics  of  the  problem  dic¬ 
tates.  Tn  either  case,  the  grid  movement  in  general  is  such  as  to  equidis- 
tribute  the  truncation  error  or,  as  is  normally  the  case,  some  measure  of  the 

L  ruination  error,  over  the  grid. 


78.  In  the  variable  number  of  grid  points  approach,  points  are  added  or 

deleted  until  some  specified  maximum  value  of  the  error  measure  allowed  over  a 

grid  interval  is  attained.  Depending  upon  the  value  specified  in  the  criterion 
to  be  satisfied,  the  variable  node  point  approach  may  yield  a  more  accurate 
solution  than  the  fixed  node  point  approach.  However,  it  will  normally  result 
in  more  difficult  coding  and  increased  computational  costs  since  the  number  of 
points  will  generally  increase.  Thus  the  most  commonly  held  view  of  what 
constitutes  adaptive  grids  is  that  a  fixed  number  of  points  are  made  to  move 

to  better  resolve  field  gradients.  A  combination  of  the  two  approaches  is, 
of  course,  also  possible. 

Equ idistr ibuted  quantities 

79.  As  in  the  fixed  number  of  grid  points  approach,  some  measure  of  the 

truncation  error  must  be  selected  to  be  equidistr ibuted  over  the  grid.  The 

grid  generator  of  Dwyer,  Smooke,  and  Kee  (1982)  is  based  upon  equidistributing 
both  the  gradient  of  the  dependent  solution  variables  given  by  as  well  as 

the  second  derivatives,  i.e.,  a  mesh  consisting  of  N  points  is  sought  such 
that 
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where  j  =  1,  2,  ...  N  -  1  and  i  =  1  ,  2  ,  ...  number  of  solution  vari¬ 
ables,  6  and  >  are  arbitrarily  selected  small  numbers,  and  max— ►/..  and 
max— ►d/.^/dx  are  determined  from  the  solution  on  the  previous  grid.  To  en¬ 
sure  adequate  smoothness  of  tin-  grid,  the  grid  spacing  is  bounded  by 


--- J  -  <  f. 
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(67) 


where  C  is  some  specified  constant  less  than  1.0  . 


Equations  65-66  in 


essence  define  the  weight  functions  that  might  be  used  in  a  variational  ap¬ 
proach  employing  a  fixed  number  of  points  or  the  |e|  in  the  attraction- 
repulsion  approach.  Both  steep  gradient  regions  as  well  as  regions  of  large 
curvature  in  any  of  the  solution  variables  would  be  resolved  by  satisfying 
Equations  65-66. 

Coniput a tional  aspects 

80.  A  new  grid  would  be  generated  by  proceeding  through  a  number  of 
steps.  The  governing  partial  differential  equations  involving  the  Z.  vari¬ 
ables  would  first  be  solved  on  an  extremely  coarse  uniformly  spaced  grid.  The 
maximum  values  of  |Z^|  an<^  jdZ^/dxj  would  then  be  computed  and  inequali¬ 
ties  such  as  those  in  Equations  65-66  evaluated.  If  either  is  not  satisfied 
in  some  interval,  a  grid  point  is  inserted  at  the  middle  of  the  interval. 

After  this  mesh  is  generated,  the  inequality  in  Equation  67  is  checked  to  en¬ 
sure  adequate  smoothness.  If  Equation  67  is  not  satisfied  in  some  interval, 
once  again  a  grid  point  is  inserted  at  the  middle  of  the  interval.  The  old 
solution  is  then  interpolated  onto  the  new  grid  and  the  steps  described  above 
are  repeated.  Dwyer,  Smooke,  and  Kee  (1982)  indicate  such  a  grid  generator 
can  be  directly  extended  to  two  dimensions  as  long  as  sharp  fronts  align  rea¬ 
sonably  well  with  one  of  the  coordinates.  If  not,  techniques  similar  to  those 
discussed  by  Anderson  and  Rai  (1982)  for  aligning  grids  along  shocks  would 
have  to  be  incorporated . 


PART  IV:  APPLICABILITY  OF  ADAPTIVE  GRIDS  IN  HYDRODYNAMIC  MODELING 


Bas i c  Concerns 

81.  The  ability  of  adaptive  grids  to  reduce  numerical  error  and  thus 
enahle  the  computation  of  more  accurate  solutions  has  been  demonstrated.  As 
previously  noted,  most  of  the  research  on  the  use  of  adaptive  grids  has  been 
in  the  aerodynamics  field  where  compressible  flow  solutions  lead  to  the  com¬ 
putation  of  large  localized  solution  jumps,  i.e.,  shock  waves.  Significant 
applications  have  also  been  made  in  flame-front  problems.  Before  considering 
particular  types  of  hydrodynamic  problems  for  which  numerical  solutions  are 
often  sought,  some  basic  concerns  with  the  application  of  adaptive  grids  in 
hydrodynamic  modeling  are  discussed. 

I n t e r p o 1 ation  prob lems 

82.  Physical  data.  Most  numerical  hydrodynamic  models  are  either  1- 
or  2-D  depending  upon  the  spatial  averaging  (Appendix  A).  Therefore  the 
geometry  of  the  water  body  is  prescribed  by  a  finite  number  of  data  points 
fixed  in  physical  space.  For  example,  in  depth-averaged  models  a  major  data 
input  is  the  water  depth,  usually  prescribed  at  the  center  of  finite  differ¬ 
ence  computational  cells.  If  the  computational  grid  is  allowed  to  move, 
obviously  a  new  water  depth  must  be  prescribed  at  the  center  of  each  new  grid 
cell.  Therefore  interpolation  of  the  original  water  depth  field  is  required 
each  time  the  grid  moves.  Although  such  an  interpolation  should  not  result  in 
additional  solution  error  since  interpolation  from  bottom  contours  is  usually 
required  to  construct  the  original  data  set,  it  will  increase  the  computational 
time.  Another  concern  along  these  lines  relates  to  bottom  roughness  coeffi¬ 
cients,  i.e.,  either  Chezy  or  Manning  coefficients.  The  normal  procedure  in 
hydrodynamic  modeling  of  free-surface  flows  is  to  calibrate  the  model  to 
known  historical  data  through  the  variation  of  primarily  the  roughness  coeffi¬ 
cient.  It  would  appear  that  such  a  calibration  would  have  to  take  place  on  a 
fixed  grid  to  generate  t  lie  appropriate  roughness  coefficient  field  for  inter¬ 
polation  to  provide  values  on  a  moving  or  adaptive  grid  in  other  applications. 

83.  Solutions  at  fixed  locations.  In  normal  free-surface  hydrodynamic 
modeling  exercises,  a  major  interest  is  in  analyzing  solution  variables,  e.g., 
the  water-surface  elevation,  as  a  function  of  time  at  some  fixed  location  such 
as  a  tide  recording  station.  If  adaptive  grids  are  used,  the  locations  of 


points  where  the  soLution  is  generated  continual Ly  move;  thus  once  again  inter¬ 
polation  at  each  time  increment  desired  for  plotting  must  be  performed  with  a 
subsequent  increase  in  the  total  cost  of  computations. 

84.  Separate  grids.  A  major  use  of  hydrodynamic  models  is  to  provide 
flow  fields  for  use  in  water  quality,  sediment  transport,  and  salinity  intru¬ 
sion  studies.  Depending  upon  the  influence  on  the  water  density  of  the  sub¬ 
stance  being  transported  and  diffused,  the  flow  and  transport  computations  may 
or  may  not  be  coupled.  In  any  case,  an  adaptive  grid  based  upon  the  equidis- 
tribution  of  some  error  measure  based  upon  the  flow  variables  may  not  be  a 
good  grid  to  minimize  numerical  errors  in  the  transport  computations  for  some 
quality  constituent.  Dwyer,  Smooke,  and  Kee  (1982)  experienced  this  problem 
in  flame  computations  involving  a  burning  spherical  particle  over  which  the 
flow  had  separated.  In  this  calculation,  both  the  flow  field  and  the  tempera¬ 
ture  gradients  had  to  be  well  resolved.  However,  the  coordinate  system  used 
in  the  temperature  computations  was  not  well-suited  for  the  flow  computations. 
As  illustrated  in  Figure  4,  Dwyer's  approach  consisted  of  using  two  different 
coordinate  systems  with  interpolation  between  the  two.  In  the  past,  hydrody¬ 
namic  modelers  have  cautioned  against  the  use  of  separate  grids  for  flow  and 
transport  computations  since  the  interpolation  of  the  flow  field,  e.g.,  total 
depth  and  velocities  in  depth-averaged  computations,  results  in  mass  not  being 
strictly  conserved.  Mass  conservation  errors  can  result  in  large  errors  in 
transport-diffusion  coinputat ions . 

85.  Another  alternative  is  to  use  one  grid  but  to  have  the  adaption  oc¬ 
cur  wherever  needed  by  any  variable.  Thus  there  would  be  regions  of  point  con¬ 
centration  due  to  the  flow  and  other  regions  due  to  the  transport.  This  would 
require  the  use  of  more  grid  points  or  else  the  mult iole-point  concentration 
regions  would  deplete  the  points  available  for  the  field  at  large.  However, 
there  would  not  be  a  need  for  interpolation  as  required  for  two  separate  grids. 

86.  Discontinuous  grid  movement.  Unless  the  time  derivatives  in  the 
governing  differential  equations  of  motion  are  transformed  and  grid  computa¬ 
tions  are  made  at  every  time-step,  interpolation  of  the  numerical  solution 
from  the  old  grid  to  the  new  grid  is  required  to  initiate  the  computations 
every  time  a  new  grid  is  computed.  Since  ttie  basic  system  of  equations  to  be 
solved  is  hyperbolic  in  character,  the  effect  of  initial  errors  dies  out 
fairly  quickly  in  stable  conputations .  Therefore  interpolation  of  an  adaptive 
grid  rn. i y  not  be  a  major  concern. 
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Computational  costs 

87.  Grid  generation.  The  typical  time  required  on  a  CRAY  computer  to 
generate  an  initial  coordinate  system  containing  about  1,000  net  points  using 

'.oiTipson '  s  coordinate  generation  code  (Thompson  et  a!.  1974),  is  approximately 
20  sec  for  100-200  iterations  of  the  Successive  Over-Relaxation  (SOR)  solution 
scheme.  This  time,  of  course,  varies  for  particular  problems  and  for  different 
degrees  of  coordinate  control.  It  should  be  remembered  that  this  generation 
code  is  based  upon  the  variational  approach  where  only  smoothness  is  involved. 
Therefore  the  computational  time  required  for  a  solution  of  the  fuLl  genera¬ 
tion  equations  given  by  Equation  45  would  be  substantially  higher.  However,  it 
also  should  be  remembered  that  the  discussion  above  relates  only  to  the  initial 
grid  generation,  where  more  than  likely  the  initial  guess  of  the  solution  is  a 
poor  one.  In  subsequent  computations  for  an  adaptive  grid,  the  initial  guess 
would  be  the  old  grid.  Normally,  there  will  not  be  major  changes  between  the 
old  and  new  grids,  so  that  perhaps  only  two  to  three  iterations  will  be  re¬ 
quired  for  each  subsequent  grid  generation.  In  other  words,  no  more  than  1  to 
2  sec  of  computational  time  wiLl  be  required  for  the  generation  of  each  new 
grid  containing  approximately  1,000  grid  points. 

88.  It  must  be  realized,  however,  that  hydrodynamic  calculations  can 
often  e:- d  ■  nd  over  long  periods  of  time.  For  example,  salinity  intrusion  com¬ 
putations  on  the  bower  Mississippi  River  by  Johnson  and  Boyd  (in  preparation) 
extended  over  4  to  6  months  with  computational  time-steps  of  15  min,  resulting 
in  perhaps  15,000  compu ta t iona l  steps  requiring  30  to  45  min  on  a  CRAY  com¬ 
puter.  If  the  grid  (x3,000  points)  had  been  computed  at  every  time-step,  an 
additional  10  to  20  hr  of  computer  time  would  have  been  required.  Obviously, 
this  could  not  be  tolerated  and  in  fact  would  not  have  been  needed.  The  move¬ 
ment  of  the  saLt  front  is  such  that  if  an  adequate  grid  had  been  used  the  grid 
would  probably  not  need  to  be  computed  more  than  once  a  day.  This  would  then 
result  in  perhaps  15  min  being  required  for  the  grid  computations,  i.e.,  per¬ 
haps  one-third  of  the  total  computational  costs  would  be  attributed  to  the  grid 
generation. 

89.  If  an  adaptive  grid  had  been  used  in  the  study  above,  the  governing 
equations,  which  were  soLved  in  a  Cartesian  form,  wou Ld  have  to  be  transformed. 
This  would  probably  result  in  an  increase  of  perhaps  two  to  three  times  the 
computational  time  cited,  i.e.,  instead  of  30  to  45  min  the  total  computer 
time  required  for  the  flow  and  salinity  computations  would  have  been  perhaps 
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1  to  2  fir.  In  this  case,  the  cost  of  generating  the  adaptive  grid  would  have 
been  perhaps  10  to  15  percent  of  the  total  cost.  For  other  more  dvnamic  free- 
surface  flow  computations,  e.g.,  storm  surge  computations,  the  percentage  of 
the  total  cost  would  be  much  higher  since  the  grid  would  need  to  be  computed 
at  virtually  every  time-step. 

90.  Time-step  restrictions.  If  ttie  time  integration  scheme  for  solving 
the  governing  partial  differential  equations  is  an  explicit  one,  the  maximum 
allowable  time-step  size  will  be  controlled  by  the  minimum  size  of  the  compu¬ 
tational  cells.  Therefore,  if  an  adaptive  grid  is  used  which  results  in  an 
extremely  smalL  grid  spacing  in  the  vicinity  of  large  gradients,  a  small,  time- 
step  will  he  required  which  in  turn  results  in  more  computer  time  for  the  flow 
computations  than  required  for  a  fixed  mesh  with  the  same  number  of  points.  0 
course  to  achieve  essentially  the  same  accuracy  with  a  fixed  grid,  a  fine  grid 
spacing  would  have  to  be  used  over  the  entire  grid.  Thus  the  small,  time-step 
would  still  he  required  hut  now  on  a  grid  containing  more  points.  Thus,  if  a 
certain  accuracy  is  required,  the  use  of  an  adaptive  grid  should  result  in 
less  computing  time.  If  implicit  time  integration  schemes  are  utilized,  there 
will  not  be  stability  restrictions  on  the  time-step. 

91.  Higher  order  versus  adaptive  grids.  Three  approaches  can  be  taken 
when  attempting  to  compute  more  accurate  numerical  solutions  of  partial  dif¬ 
ferential  equations.  Obviously,  the*  number  of  grid  points  can  be  increased. 
However,  if  the  number  ot  grid  points  remains  fixed  either  higher  order  solu¬ 
tion  schemes  or  adaptive  grills  with  lower  order  schemes  must  be  employed.  A 
good  exanpie  of  a  higher  order  scheme  that  filters  out  dispersive  effects  is 
the  f  lux-correct  *>il  transport  scheme  of  ZaLesak  (  1979). 

'•>2.  In  its  simplest  terms,  flux-corrected  transport  constructs  the  net 
transportive  flux  point  by  point  as  a  weighted  average  of  a  flux  computed  by  a 
low-order  scheme  and  a  fmx  computed  by  a  high-order  scheme.  The  crucial  step 
is  the  determination  of  the  weighting  which  is  related  to  the  flux  limiting 
step.  Basically,  an t id i f f us ive  fluxes  are  computed  as  the  difference  between 
the  high-order  flux  and  the  Low-order  flux  and  are  then  limited  in  such  a  way 
that  the  fieid  is  free  ot  artificial  extrema.  An  example  of  the  superiority 
of  the  scheme  over  a  iow-order  scheme  is  presented  in  Figure  5. 

93.  It  is  rather  lifticult  to  compare  accuracy  and  computer  time  re¬ 
quirements  of  higher  order  schemes  versus  the  use  of  adaptive  grids  with  low- 
e  ha  pes  without  <■,),!  iig,  both  wit' in  the  same  computer  code  and  applying 


each  to  the  same  problem.  However,  experiences  cited  in  the  Literature  for 
each  approach  have  been  noted.  Dortch*  found  that  the  use  of  ZaLesak's  flux- 
corrected  transport  scheme  increased  the  computer  time  by  a  factor  of  2.7  over 


the  use  of  a  donor  ceLL  scheme  in  a  2-')  LateraLLy  averaged  dye  concentration 
computation.  Approximately  four  times  as  many  ceLLs  were  required  in  the 
donor  cell  computation  to  achieve  the  same  accuracy  as  those  obtained  with  the 
fin;:  corrected  transport  computations.  Braokb i I i ' s  2-D  computations  of  the 
steady  advective  diffusion  equation  showed  that  approximately  10  times  as 
many  cells  in  a  fixed  grid  computation  were  required  to  achieve  the  same  ac¬ 
curacy  as  that  obtained  with  an  adaptive  grid.  No  comparison  of  computer  time 
was  given. 

Hydrodynamic  Applications 

94.  In  the  discussion  below,  several  types  of  problems  routinely  ad¬ 
dressed  by  numerical  hydrodynamic  models  are  noted.  In  some  cases,  adaptive 
grids  wouLd  be  of  Little  use;  whereas  in  others,  a  significant  improvement  in 
modeling  capability  could  be  realized  through  the  use  of  adaptive  grid  tech¬ 
niques.  Of  the  grid  generation  techniques  developed  to  date,  it  is  believed 
that  the  variational  approach  of  BrackbiLl  (1982)  offers  the  most  promise; 
therefore  some  appropriate  weighting  functions  for  such  a  grid  generator  are 
noted . 

Flood  routing 

l>5.  The  computation  of  l  loods  on  natural  rivers  involves  unsteady  1-D 
computations  for  the  water  surface  and  the  flow  rate,  i.e.  discharge.  Although 
a  wave  i -  computed,  it  is  a  long-period  wave  (Figure  6)  with  a  wavelength 
norma i » v  on  the  order  of  hundreds  of  miles.  Therefore,  in  general,  many  grid 
points  nro  contained  within  a  wavelength  and  numerical  error  is  small.  For 
siK'h  pro'M"'s,  no.  adaptive  grid  would  he  of  little  use. 

Dam  breaks 

9b.  (Ver  the  past  five  years,  a  major  modeling  activity  within  the  Corps 
has  been  the  computation  of  flooding  resulting  from  hypothetical  failures  of 
Corps-regulated  reservoirs.  In  such  computations,  there  is  basically  a  large 
single  wave  (Figure  7)  moving  down  the  valley  below  the  dam.  Major  interest 


*  Personal  communication,  M.  Dortch,  1984,  US  Army  Engineer  Waterways  Ex¬ 
periment  Station,  Vicksburg,,  Miss. 


is  in  predicting  arrival  times  as  well  as  maximum  water-surface  elevations 
at  various  locations.  Theoretically,  the  use  of  a  1-D  adaptive  grid  in  such 
computations  should  esult  in  increased  accuracy  for  a  fixed  number  of  points. 
Of  course  the  problem  of  interpolation  of  channel  cross-sectional  geometry  data 
would  add  to  the  computational  costs.  In  addition,  if  the  cross  section  is 
broken  into  a  conveyance  plus  a  storage  section,  the  storage  volume  on  the 
overbank  is  the  grid  spacing  times  the  storage  cross-sectional  area.  As  the 
dam-break  peak  moves  downstream,  the  grid  upstream  will  become  sparse.  As  a 
result,  it  is  easy  to  visualize  continuity  problems  related  to  the  storage 
volume  occurring  in  the  computations.  Therefore,  in  a  practical  sense,  it  may 
be  better  to  make  computations  on  a  grid  containing  enough  points  to  yield  the 
d  e s i r e d  a  c  c ur a e y . 

97.  If  an  adaptive  grid  is  employed  in  a  dam-break  model,  a  weight 
function  based  upon  the  water-surface  elevation  is  suggested.  Using  Equa¬ 
tion  61,  the  suggested  weight  function  becomes 

w  =  (1  +  i;2  K  )  yj  1  +  a2h2  (68) 


1  +  h 


where 

h  =  water-surface  elevation 

i,  8  =  positive  parameters  determined  from  numerical  experimentation 
S_u rg_e  c o mputat  ions 

98.  One-dimensional .  The  sudden  release  of  large  volumes  of  water 
either  from  hydropower  peaking  operations  or  perhaps  as  a  result  of  lock  opera¬ 
tions  can  create  monoclinical  gravity  waves  that  are  normally  referred  to  as 
surges.  Positive  surges  that  occur  downstream  of  such  operations  elevate  the 
water  surface,  whereas  negative  surges  that  occur  upstream  depress  the  water 
love  1 .  The  study  of  surges  is  important  because  of  their  possible  adverse 
influence  upon  navigation,  as  well  as  the  sudden  large  loads  they  can  impose 
upon  lock  miter  gates  and  appurtenances.  Surges  of  large  amplitude  cannot  be 
tolerated  in  practical  navigation  channels;  therefore  accurate  predictive  com¬ 
putations  are  essential. 


‘,l’  •  111*-'  computation  ot  such  stirpes  involves  solving  the  1-1)  unsteady 

flow  eiiu. .t  ions  presented  in  Appendix  A.  Many  such  models  have  been  developed 
(Johnson,  in  preparation).  However,  a  deficiency  of  all  such  models  is  the  nu- 
■ :  i>  i  i  >  .p<  ion  and  dissipa  tion  that  is  inherent  when  attempting  to  compute 

waves  that  are  resolvcii  over  on  1 \  a  lew  grid  points,  e.g.,  < 1 0  .  Such  problems 
arc  ideal  lv  suited  tor  adaptive  griddine,  especially  if  a  single  surge  is 
created,  however,  hviiropeaking  operat  ions  often  it. ate  more  than  one  surge, 
figure  8  is  an  example  of  typie.al  surges  created  downstream  of  Barkley  Dam. 
Althou.Ji  multiple  surges  will  not  introduce  p..rt  ieuiar  problems  with  applying 
adapt i ve  ..rids,  more  total  points  will  he  required  so  that  there  will  he  enough 
points  to  get  eoneent  rat  ion  in  each  surge  without  depleting  the  field  else¬ 
where.  As  a  surge  passes  out  the  downstream  boundary,  of  course  more  points 
become  available  elsewhere.  The  more  surges  expected  in  the  field  at  any  one 
time,  the  greater  will  he  the  eomput at i onal  time  since  the  total  number  of 
points  will  be  dependent  upon  the  number  of  surges. 

100.  As  with  t hr*  dam-break  problem,  interpolation  of  the  cross-sectional 
geometry  tables  would  be.  required.  However,  the  problems  connected  with  over¬ 
bank  storage  would  probably  not  be  as  prevalent  since  overbank  areas  would  not 
be  as  likely  to  be  flooded.  Since  water  velocity  and  surge  height  are  major 
considerations  in  assessing  the  impact  of  surges  on  navigation,  adaption  could 
probably  he  based  upon  either.  Once  again  the  weight  function  given  by  Equa¬ 
tion  b8  is  suggested  if  adaption  is  based  upon  the  surge  height;  whereas  the 
wat er-sur face  elevation  would  he  replaced  in  Equation  <18  hv  the  average  cross- 
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Short  waves  in  harbors 

102.  The  behavior  of  short  waves  in  shallow  water  is  a  major  concern  of 
coastal  engineers.  The  design  of  harbors  and  other  facilities  necessitates  a 
detailed  knowledge  of  short-wave  characteristics  since  these  waves  attack 
harbor  structures  while  being  refracted,  diffracted,  and  reflected  through  the 
harbor.  An  interesting  computation  by  Brackbill  (1982)  involving  aerodynamic 
wave  reflections  on  an  adaptive  grid  is  presented  in  Figure  10.  Such  waves 
also  influence  navigation  in  the  harbor  area  directly  and  through  the  currents 
tiiev  can  induce. 

101.  In  order  to  compute  the  behavior  of  such  waves,  additional  terms 
must  lie  added  to  the  depth-averaged  equations  presented  in  Appendix  A  to  ac¬ 
count  for  the  effect  of  vertical  accelerations.  The  introduction  of  these 
terms  presents  problems  in  that  solutions  of  the  equations  are  quite  sensitive 
to  t he  numerical  errors  introduced  by  standard  numerical  techniques  unless  per¬ 
haps  10  to  20  grid  [joints  are  used  per  wavelength.  However,  from  a  practical 
standpoint,  there  is  a  need  to  describe  any  short  wave  with  only  a  few  grid 
points.  As  with  the  surge  computations,  it  would  appear  that  adaptive  grids 
would  be  the  solution  to  the  problem.  However,  as  illustrated  in  Figure  11, 
many  short  waves  are  genera  1 1 y  contained  within  the  computational  grid.  Thus 
it  is  doubtful  that  an  adaptive  grid  would  vield  much  benefit  since  the  total 
number  el  points  required  for  good  concentration  over  each  wave  would  probably 
be  about  the  same  as  that  required  with  a  fixed  grid.  A  better  approach  might 
he  to  employ  some  higher  order  scheme,  such  as  the  flux-corrected  transport 
previously  discussed,  applied  on  a  grid  containing  as  many  points  as  economi¬ 
cal  I v  possible. 

I  i  i.l a  1  c  I  rc:u  I  a  t  i  ons 

I0'i.  Vertically  ave raped  numerical  models  are  often  emp loved  to  compute 
t  low  t  ields  in  relatively  shallow,  we  11 -nixed  estuaries  and  coastal  areas  for 
use  in  ■■.siler  quality  and  sedinieiitut  ion  studies.  Examples  are  the  models  hv 
[’■uf  let  f  lb, sill  and  Morton  and  King  (  1  ‘>77 )  .  Although  the  propagation  of  a  tidal 
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part  it'll  lnr  problems  sinrr  several  grid  points  lie  witliin  a  wavelength.  Al¬ 
though  large  gradients  in  the  (low  do  oeeur  when  moving  from  deep  areas  to 
relatively  shallow  ones,  e.g.,  from  navigation  channels  to  the  shallow  over- 
ohs  i,l  i  a.  re  12),  it  is  not  believed  that  in  general  a  moving  grid  would  he 
ot  use  since  basieallv  the  regions  of  high  gradients  are  known  and  don't 
arbitrarilv  develop  in  the  computational  grid.  However,  it  is  believed  that 
adaptive  grid  techniques  would  be  of  great  use  in  numerical  1 v  generating  the 
initial  grid. 

103.  Navi  station  channels.  In  Thompson's  grid  generation  code,  control 
(unctions  can  he  prescribed  to  force  a  small  grid  spacing  in  particular  re¬ 
gions.  However,  the  control  is  relative  to  either  f  or  n  lines  in  the  com 
putational  field.  I'nless  the  line  to  which  attraction  is  desired  lies  on  a 
boundary,  its  position  in  physical  space  is  not  known  until  after  the  grid  has 
been  generated.  Attraction  to  specific  lines  in  space  is  more  difficult  and 
more  expensive  and  has  not  proved  to  he  as  successful.  In  other  words,  cur¬ 
rently  it  is  no t  possible  with  this  code  to  force  a  fine  grid  spacing  along 
winding  navigation  channels  fixed  in  phvsical  space.  The  generation  of  an 
adaptive  grid  which  would  only  he  generated  once  should  provide  a  solution  to 
this  problem.  Equation  A3  would  constitute  the  grid  generator  with  the 
weight  function,  w  ,  being  prescribed  so  as  to  force  a  small  grid  spacing 
where  w  is  large  and  a  large  spacing  where  w  is  small.  Since  navigation 
channels  are  defined  by  their  increased  depth,  the  appropriate  weight  function 
should  he  the  water  depth  itself,  i.e., 

w ( x , y )  —  h ( x , v )  (70) 

where  h(x,v)  is  the  total  water  depth  relative  to  some  mean  water  level. 

This  would  have  the  effect  of  generating  a  grid  with  the  same  volume  in  each 
vertical  column. 

10b.  Flooding  and  drying.  As  computations  are  made  over  a  tidal  cycle, 

quite  often  there  is  a  need  to  allow  for  flooding  of  tidal  flats  on  the  flood 

cycle  of  the  tide,  with  a  subsequent  drying  of  those  areas  on  the  ebb  por¬ 
tion.  In  other  words,  the  horizontal  water  boundary  moves  over  a  tidal  cycle. 

One  could  think  in  terms  of  using  adaptive  gridding  to  handle  such  moving 
boundaries.  An  example*  of  a  moving  boundary  adaptive*  grid  is  presented  in 
Figure*  13.  Hi  iwe*vo  r ,  if  the  am  a.  flooded  and  dried  are  essentially  internal 


areas,  e.g.,  an  island  that  appears  and  then  disappears,  some  technique  for 
allowing  points  to  collapse  into  a  single  point  when  the  area  is  flooded  would 
he  required.  In  addition,  moving  points  around  to  better  represent  the  water 
boundary  would  probably  be  to  the  detriment  of  the  channel  resolution.  As  a 
result  of  these  anticipated  problems,  it  is  believed  that  the  best  approach 
would  be  to  use  adaptive  grid  techniques  to  generate  the  initial  grid,  with 
grid  lines  generally  following  bottom  elevation  contours.  Wetting  and  drying 
could  then  be  carried  out  on  this  grid  as  is  currently  done  on  Cartesian  grids, 
i.e.,  bv  removing  or  adding  cells  from  the  computations  by  checking  water- 
surface  levels. 

107.  Finite  element  stability.  On  occasion,  numerical  experience 
confirms  that  the  Norton-Ring  (1977)  finite  element  model  used  within  the 
Hydraulics  Laboratory  can  become  unstable  midway  in  a  model  execution  due  to 
the  particular  grid  used.  The  current  "cure"  is  to  redistribute  nodes  or 
channel  topography,  in  a  manner  based  on  experience,  and  execute  the  code 
again.  For  a  mathematically  consistent  algorithm,  the  approximation  error  is 
always  bounded  by  some  measure  (size)  of  the  mesh  (grid),  i.e.,  the  error  is 
controlled  by  refining  (perhaps  locally)  the  grid.  An  improvement  to  the 
utilization  of  the  codes,  regarding  tidal  cycle  restarts,  could  result  from 
extension  of  the  theory  to  include  node  motion.  An  energy  norm  distribution, 
say  in  h  ,  the  free-surface  height,  could  detect  regions  of  incipient  non¬ 
smoothness  (growing  instability)  and  hence  pinpoint  regions  in  which  a  refined 
grid  could  be  appropriate.  Within  the  code  itself,  this  amounts  to  adding 

the  grid  convection  velocity  terms,  :X}  in  Equation  33. 

St  rat  L  f led _ flow 

108.  The  computation  of  flow  fields  in  bodies  of  water  where  the  den¬ 
sity  varies  over  the  depth  often  involves  internal  density  currents.  Such 
stratified  flow  refers  to  motions  involving  fluid  masses  of  the  same  phase.  A 
heavier  liquid  flowing  beneath  a  lighter  liquid,  or  a  heavier  gas  moving  under 
a  lighter  gas,  will,  he  subject  to  gravitational  effects  that  depend  upon  the 
differences  between  the  two  densities.  Figure  14  illustrates  such  a  flow 
field  for  a  cold-water  inflow  near  the  bottom  of  a  flume.  The  internal  motions 
in  reservoirs  due  to  temperature  variations  and  salinitv  intrusion  into  rela¬ 
tively  deep  estuaries  are  examples  of  commonly  computed  stratified  flow  fields. 

109.  Ihe  salinity  intrusion  study  bv  Johnson  and  Boyd  (in  preparation) 
on  the  Lower  Mississippi  River  is  a  good  example  of  numerical  computations 


involving  density  otto  ts.  !  he  cavern i ng  equations  solved  are  the  laterally 
averaged  equations  lor  both  the  i low  t ield  and  the  salinity  field  presented  in 
Appendix  A.  Since  the  salinity  computations  involve  the  computation  of  a  mov- 
'  ng  salt  wedge  with  a  sharp  inierlact  separating  the  sal t  water  from  the  fresh 
water  flowing  down  the  river  above  the  wedge,  adaptive  gridding  should  enable 
the  computation  of  a  much  sharper  front. 

110.  As  in  the  Dwyer,  Smooke,  and  Kee  (1982)  flame  problem,  it  may  be 
that  separate  grids  would  have  to  he  employed  for  the  flow  and  salt  computa¬ 
tions.  However,  due  to  the  interpolation  problems  previously  discussed,  be¬ 
fore  more  than  one  grid  is  employed,  the  use  of  a  single  grid  should  first  be 
attempted.  Since  gradients  in  both  the  horizontal  velocity  and  the  salt  con¬ 
centration  are  large  in  the  vicinity  of  the  saltwater-freshwater  interface,  an 
appropriate  weighting  function  to  generate  a  single  grid  should  involve  both. 
For  example,  the  following  weight  function  might  be  applied: 

W ( x , z )  =  rtyj  1  +  VtT  +  (1  -  a)  yj  1  +  Vs2  (71) 

where 

a  =  constant  with  values  between  0  and  1 
U  =  horizontal  velocity 
s  =  salinity 
H y d r a u 1 i c  jumps 

111.  When  bottom  slopes  change  from  being  steep  enough  to  support  super¬ 
critical  flow  to  milder  slopes,  subcritical  flow  conditions  are  arrived  at 
through  the  formation  of  a  hydraulic  jump.  This  is  completely  analogous  to 
the  formation  of  shock  waves  in  compressible  aerodynamic  flow  fields.  Such 
flow  conditions  are  rare  in  nature  but  can  occur  downstream  of  spillways  at 
flow-control  structures.  The  analysis  of  this  free-surface  hydrodynamic 
problem  is  normally  accomplished  with  small-scale  physical  models.  However, 
since  the  formation  of  shock  waves  is  now  routinely  computed  in  aerodynamics, 
the  numerical  computation  of  the  formation  of  hydraulic  jumps  should  also  be 
possible.  In  fact.  Figure  15  presents  the  results  of  a  preliminary  computa¬ 
tion  of  a  hydraulic  jump  by  Baker*  using  a  finite  element  code.  Obviously, 

*  Personal  communication,  A.  .1.  Baker,  1983,  University  of  Tennessee, 
Knoxville,  Tenn. 
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use  of  an  adaptive  grid  in  such  computations  would  be  of  great  benefit;  and  the 
experience  of  compressible  aerodynamic  modelers  would  be  directly  applicable. 

In  a  similar  fashion  to  the  pressure  gradient  being  used  in  the  weighting  func¬ 
tion  in  many  aerodynamic  adaptive  grid  applications,  the  appropriate  weighting 
function  should  probably  be  based  on  the  gradient  of  the  water  depth. 


PART  V :  SUMMARY 


112.  Since  numerical  solutions  of  the  partial  differential  equations 
governing  the  hydrodynamics  of  bodies  of  water  are  obtained  on  grids  that 
represent  a  discretization  of  the  continuous  physical  domain,  errors  that  are 
a  function  of  t lie  discretization  invariably  arise.  Historically,  such  errors 
have  been  controlled  in  one  of  two  ways.  Either  a  complex  higher  order  solu¬ 
tion  scheme  is  employed,  or  an  extremely  fine  resolution  grid  is  constructed 
over  the  physical  domain  so  that  simple  low-order  schemes  can  be  used.  From  a 
programming  and  computational  efficiency  viewpoint,  low-order  schemes  are  more 
desirable.  However,  increasing  the  number  of  computational  points  is  often  an 
inefficient  approach  since  portions  of  the  domain  may  contain  small  gradients 
of  the  solution  variables  and  a  fine  grid  is  not  required  there.  The  solution 
to  this  problem  is  to  employ  a  grid  that  is  free  to  move  as  the  solution  is 
computed  such  that  a  fine  grid  spacing  occurs  in  regions  of  large  solution 
gradients  with  a  corresponding  larger  grid  spacing  in  small  gradient  regions. 
Such  grids  are  referred  to  as  adaptive  grids. 

Cone lusions 


113.  Adaptive  grids  can  be  implemented  in  both  finite  element  solution 
schemes  and  finite  difference  solution  schemes  applied  on  boundary-fitted  co¬ 
ordinates.  The  most  natural  implementation  involves  transforming  the  time 
derivative  in  the  governing  flow  equations  so  that  computations  are  performed 
in  a  fixed  transformed  grid  with  the  influence  of  the  grid  movement  in  the 
physical  plane  entering  the  computations  through  terms  involving  the  speed  of 
the  grid  points. 

114.  The  usual  definition  of  an  adaptive  grid  is  one  containing  a  fixed 
number  of  points,  with  the  movement  of  those  points  dictated  by  the  physics  of 
the  problem  being  solved.  Several  methods  for  generating  such  grids  exist  and 
there  are  advantages  and  disadvantages  to  each.  However,  all  basically  are 
concerned  with  the  equ id ist r ibu t ion  of  some  measure  of  the  solution  error  over 
the  grid.  The  approach  that  seems  to  be  the  most  logical,  and  that  offers  the 
greatest  inherent  control  over  the  grid,  is  a  variational  approach  developed 

by  Brackbill  (1982).  A  functional  whose  value  is  given  by  an  integral  equation 
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involving  integrals  related  to  grid  smoothness  and  orthogonality,  as  well  as 
the  concentration  of  grid  points,  is  minimized  to  yield  differential  equations 
for  computing  the  grid  point  locations.  Therefore  different  emphasis  can  be 
placed  upon  smoothness,  orthogonality,  and/or  point  concentration  merely  by 
changing  the  value  of  three  constants.  A  major  question  is  what  to  use  for 
the  weight  function  that  controls  the  grid  adaption. 

115.  Several  hydrodynamic  problems  that  are  often  addressed  with  numer¬ 
ical  models  have  been  discussed  with  respect  to  the  usefulness  of  obtaining 
solutions  on  adaptive  grids,  and  suggestions  for  weight  functions  have  been 
given.  In  general,  if  the  flow  solution  contains  perturbations  that  either 
develop  and/or  move  in  arbitrary  regions  of  the  physical  domain,  adaptive  grid- 
ding  will  be  of  greatest  benefit,  e.g.,  the  computation  of  surges  as  a  result 
of  either  locking  or  hydropower  peaking  operations,  and  the  computation  of  hy¬ 
draulic  jumps  and  concentration  fronts.  However,  if  many  perturbations  exist 
in  the  domain  at  any  particular  time,  adaptive  gridding  probably  offers  no 
advantage.  The  propagation  of  short  waves  in  coastal  areas  is  an  example. 

Re c  ommenda t ions 

116.  Much  of  the  numerical  modeling  of  free-surface  hydrodynamics 
within  the  Corps  of  Engineers  involves  the  depth-averaged  computation  of  cir¬ 
culation  patterns  in  estuaries  and  coastal  areas.  In  such  problems,  a  major 
concern  is  the  representation  of  winding  navigation  channels.  Although  time- 
dependent  adaptive  grids  don't  in  general  appear  to  be  required  in  the  solu¬ 
tion  of  such  problems,  a  grid  generator  based  upon  the  Brackbill  variational 
approach  should  be  an  ideal  generator  for  the  initial  grid.  Resolution  in  the 
navigation  channel  could  be  achieved  by  adapting  to  the  water  depth  with  a 
resulting  larger  spacing  in  the  nonchannel  areas.  Such  a  grid-generating  model 
is  required  before  the  finite  difference  models  of  Johnson  (1980)  and  Sheng 
and  Hirsh  (1984),  based  upon  boundary-fitted  coordinates,  can  be  generally  ap¬ 
plied.  In  addition,  it  should  also  make  the  construction  of  finite  element 
grids  a  much  more  automated  process.  Therefore  it  is  recommended  that  a  re¬ 
search  effort  be  initiated  to  develop  an  estuarine/coastal  fixed-grid  generator 
model  based  upon  adaptive  grid  techniques. 

117.  Although  a  general  adaptive  grid  is  not  needed  in  tidal  circula¬ 
tion  models,  numerical  experience  confirms  that  the  depth-averaged  finite 


rlomeiu  model  (RMA-2)  often  employed  by  the  Corps  can  become  unstable  midway 
in  a  model  execution  it  a  grid  that  does  not  adequately  resolve  flow  gradients 
in  localized  regions  is  employed.  Research  into  how  to  best  incorporate 
.d apt  ive  .'rid  techniques  into  the  existing  model  to  automatically  make  minor 
local  changes  in  the  finite  element  grid  will  provide  the  potential  for  tre¬ 
mendous  savings  in  computer  time  and,  therefore,  is  highly  recommended. 

118.  At  the  present  time,  small-scale  physical  models  are  employed  to 
address  the  design  of  spillways  where  hydraulic  jumps  can  occur.  The  numerical 
techniques  employed  in  the  development  of  models  for  computing  the  transition 
flow  from  supersonic  to  subsonic  (and  shock  waves  that  occur)  should  be  uti¬ 
lized  to  develop  hydrodynamic  transition  flow  models.  Adaptive  grids  will 
enable  such  computations  to  be  made  more  accurately  at  less  cost. 
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Figure  1.  Transformation  of 
physical  domain  to  computa¬ 
tional  (note  treatment  of 
is  1  and) 


Figure  2.  Example  of  an  1-D  adap 
tive  grid  (from  White  1982) 
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a.  Coordinate  system  for 
temperature  computations 
and  isotherm  distribution 


b.  Coordinate  system  for 
flow  computation  and 
vorticity  distribution 


Figure  4.  Example  of  multiple  coordinate  systems 
(from  Dwyer,  Smooke ,  and  Kee  1982) 
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Figure  ‘3.  Idealized  comparison  of  tl ux-correc ted 
transport  and  donor  cell,  (from  Johnson,  Thompson, 
and  Stein  1482') 
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igure  7.  Outflow  hydrograph  from  Teton  Dam 
failure  (from  Fread  1978) 
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b.  Downstream  surge 


Figure  8.  Surge  resulting  from  hydropower  peaking 
at  Barkley  Dam  (from  Johnson  1978) 


Figure  9.  Observed  (squares)  and  computed  (solid)  water 
levels  for  Hurricane  CamiLle  (from  Wanstrath  1977) 


Figure  12.  Current  pattern  in  Coos  Bay 
(from  Butler  1978) 


Figure  15.  Hydraulic  jump  computations  (from  Baker  1983) 


APPENDIX  A:  HYDRODYNAMIC  EQUATIONS  A..U  APPROXIMATIONS 


Approximations 


1.  The  Navier-Stokes  equations  express  the  conservation  of  mass  and 
momentum  of  a  flow  field  and  are  the  basic  governing  equations  for  the  solution 
of  any  fluid  dynamics  problem.  These  equations  written  in  tensor  notation  are 


Continuitv : 
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c)t  CX. 
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=  0 
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Momentum : 
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where 


P  =  water  density 
t  =  t ime 


u^jU^jU^  =  tensor  notation  for  velocity 


x.,x.  =  tensor  notation  of  spatial  coordinate 
i  J 


P  =  pressure 


=  acceleration  due  to  gravity 


=  cyclic  tensor 

ijk 


Slj  =  Coriolis  parameter 

T..  =  laminar  stress  tensor 
ij 

In  addition  to  the  above  equations,  a  conservation  of  energy  equation  must 


also  be  written  for  fluid  flow  problems  with  thermal  effects.  With  the  as¬ 
sumption  of  a  constant  specific  heat  and  with  the  neglection  of  viscous  dis¬ 
sipative  effects,  the  energy  equation  can  be  written  as  the  following  transport 
equation  for  temperature  T  : 


Energy : 


IT 

at 


3 (Tu.) 
3x  . 


I 


sinks 


(A3) 


where  D. .  is  the  diffusivity  coefficient.  This  equation  states  that  the  tem- 
i.l 


perature  can  change  as  a  result  of  advection  by  the  flow  field,  the  molecular 


diffusion,  and  the  actions  of  any  sources  and  sinks  of  heat.  As  a  matter  of 
fact,  this  same  equation  applies  to  the  transport  of  any  constituent  0  , 
where  0  replaces  T  in  the  equation  and  appropriate  sources,  sinks,  and 
boundary  conditions  are  prescribed.  For  example,  in  the  numerical  modeling  of 
the  hydrodynamics  of  an  estuary,  a  similar  transport  equation  for  the  salinity 
would  be  required. 

2.  One  additional  equation  remains  to  be  written  in  order  to  close  the 
system.  An  equation  of  state  expressing  the  density  as  a  function  of  the 
temperature  and  pressure  (and  salinity  in  estuarine  modeling)  must  be  employed 

Equation  of  State:  p=p(T,P)  (A4) 

3.  With  the  closure  of  the  system,  six  equations  exist  to  be  solved  for 
the  six  unknowns — density  p  ;  three  velocity  components  u  ,  v  ,  w  ;  pres¬ 
sure  P  ;  and  temperature  T  . 

4.  The  above  equations  written  with  molecular  values  of  viscosity  and 
diffusivity  are  only  applicable  to  laminar  flows.  Following  Reynolds,  the 
approach  taken  to  make  the  equations  applicable  to  turbulent  flows  is  to  as¬ 
sume  that  each  dependent  variable  is  composed  of  an  average  time-varying 
component  plus  a  small  randomly  varying  component  abort  the  average  value. 
Integration  over  an  increment  of  time  produces  a  set  of  equations  of  the  same 
form;  however,  these  equations  are  now  written  with  the  time-averaged  com¬ 
ponents  as  the  dependent  variables  plus  additional  terms  involving  products  of 
the  turbulent  fluctuations. 

5.  An  important  problem  in  numerical  hydrodynamic  models  is  the  treat¬ 
ment  of  these  turbulence  terms.  Most  existing  models  utilize  the  concept  of 
an  eddy  viscosity.  Eddy  viscosity  models  are  relatively  easy  to  use  and  can 
often  give  reasonable  results  if  sufficient  data  are  available  to  establish 
the  validity  of  the  model  parameters.  However,  under  flow  conditions  where 
the  turbulence  is  generally  not  in  equilibrium  with  the  mean  flow  gradients, 
e.g.,  highly  oscillatory  flows;  models  that  resolve  the  time  rate  of  change, 
convection,  diffusion,  production,  and  dissipation  of  turbulence  are  required. 
These  models  are  known  as  Reynolds  stress  models  and  involve  the  solution  of 
dynamic  equations  for  the  second-order  turbulent  correlations  along  witli  the 
equations  for  the  mean  flow  variables.  With  this  type  of  model,  a  universal 
set  of  constants  for  a  wide  variety  of  flow  conditions  can  be  used. 


6.  Subject  to  the  assumptions  made  in  their  derivation,  the  time- 
averaged  equations  are  applicable  to  any  turbulent  fluid  dynamics  problem. 

An  approximation  usually  made  when  applying  the  equations  to  hydrodynamic 
problems  was  first  pointed  out  by  Boussinesq.  When  variations  of  temperature 
are  small  (At  <_  +10°C) ,  variations  in  density  will  be  less  than  one  percent. 
These  small  variations  can  be  ignored  in  general  with  one  exception.  The 
variability  of  density  in  the  gravitaional  term  cannot  be  ignored.  Hence,  p 
is  treated  as  a  constant  in  all  places  except  the  body  force  term. 

7.  An  assumption  that  is  usually  made  in  hydrodynamic  models  is  that 
vertical  accelerations  are  negligible  when  compared  with  the  gravitational  ac¬ 
celeration.  Neglecting  viscous  terms  also,  the  vertical  momentum  equation 
reduces  to 

I  =  <A5> 

where  z  is  the  vertical  spatial  coordinate  and  is  positive  downward.  Equa¬ 
tion  A5  states  that  the  pressure  is  hydrostatic. 

8.  When  considering  the  coupling  of  the  thermodynamics  and  the  hydro¬ 
dynamics  of  a  water  body,  a  disuinction  must  be  made  between  convective  and 
"auasi-static"  models.  Convective  models  retain  the  complete  vertical  momentum 
equation  and  can  simulate  in  full  detail  the  convective  overturning  of  unstable 
water  masses,  such  as  the  upwelling  of  cells  of  warm  water  or  perhaps  the 
plunging  of  a  cold-water  inflow.  Quasi-static  models  eliminate  vertical  ac¬ 
celerations  due  to  buoyancy  effects,  which  precludes  the  explicit  treatment  of 
free  convection  associated  with  unstable  stratification.  Convective  over¬ 
turning  can  only  be  handled  as  mixing  along  the  vertical  axis.  Generally, 
models  with  the  hydrostatic  pressure  assumption  require  far  less  computer  time 
than  the  fully  convective  models. 

9.  A  spatial-averaging  technique  similar  to  the  turbulent  time-averaging 
technique  is  used  to  yield  either  1-  or  2-1)  model  (Ward  and  Espey  1971).*  The 
spatial  integration  yields  a  set  of  equations  with  the  time-averaged  and  spa¬ 
tially  averaged  components  of  the  flow  and  thermodynamic  variables  as  dependent 
variables  plus  terms  involving  products  of  the  spatial  fluctuations.  As  usu¬ 
ally  done  with  the  turbulent  fluctuations,  eddy  coefficients,  referred  to  as 


*  References  are  at  end  of  main  text 


eddy  dispersion  coefficients  to  distinguish  them  from  the  turbulent  eddy 
diffusion  coefficients,  are  employed  to  close  the  governing  equations. 

10.  Many  depth-averaged  models  have  been  developed  over  the  past 

20  years  or  so;  whereas  width-averaged  models  have  only  been  developed  over 
the  past  8  to  10  years.  If  the  spatial  integration  is  performed  over  the  com¬ 
plete  cross  section,  a  1-D  model  with  variations  allowed  only  in  the  longi¬ 
tudinal  direction  results.  The  governing  sets  of  equations  for  both  1-  and 
2-1)  free-surface  hydrodynamic  models  are  given  below. 

One-Dimensional  Cross-Sectionally  Averaged 

11.  For  the  case  of  1-D  open-channel  flows  within  rigid  boundaries,  the 
flow  behavior  can  be  adequately  described  by  the  Saint-Venant  partial  differ¬ 
ential  equations  of  unsteady  flow.  These  equations,  which  are  presented  be¬ 
low,  can  be  derived  by  integrating  Equations  A1  and  A2  over  a  cross  section. 

Flow  Continuity:  ~  ~  =  q„  (A6) 

dx  ot  i 

Momentum:  +  V  ~  (6Q)  +  3V  ^  -  6V2B 

Ot  dx  Ox  Ox 

+  gA  =  gA  (Sx  -  Sf)  +  BvV  (A7) 

12.  Equations  A6  and  A7  comprise  the  set  of  equations  governing  the 
motion  of  water  in  open  channels  in  a  1-D  sense  and  involve  three  unknowns, 
namely,  the  flow  discharge  Q  ,  the  flow  depth  y  ,  and  the  frictional  slope 
Sj.  .  Other  variables  such  as  the  lateral  inflows  and  geometry  data  are  ex¬ 
pected  to  be  known.  To  achieve  closure  of  the  system,  an  additional  relation 
is  required.  This  is  provided  by  Manning's  equation  which  relates  the  fric¬ 
tion  slope  to  the  flow  and  channel  characteristics 


s  =  -i  - oLoL 

f  2.25A2R4/,:3 


(A8) 


With  this  additional  relation,  one  can  then  solve  for  the  basic  unknowns  Q 
and  v  .  Variables  in  the  equations  are  defined  as: 


A  = 


x 


B  = 

A  = 

K  = 
e 

n  = 

Q  = 

q?  = 
R  = 

S  = 
x 


t  = 
V  = 
x  = 

y  = 

z  = 

e  = 
= 

3/Dx  = 


total  cross-sectional  area  of  channel 

derivative  of  A  with  respect  to  channel  distance  at  a  constant 
flow  depth 

top  width  of  water  surface 
acceleration  due  to  gravity 
coefficient  in  eddy  head  loss  term 
coefficient  in  Manning's  equation 
flow  discharge,  in  cfs 
lateral  inflow  rate  of  water 
hydraulic  radius 
slope  of  channel  bed 
friction  slope 
t  ime 

average  flow  velocity 

distance  along  the  channel 

depth  of  water  in  the  channel 

elevation  of  channel  bed 

momentum  correction  factor 

derivative  with  respect  to  time 

derivative  with  respect  to  channel  distance 


Two-Dimensional  Vertically  Averaged 

13.  An  integration  of  the  time-averaged  3-D  equations,  with  T  replaced 
by  the  salinity  s  ,  yields  a  set  of  equations  with  the  time-averaged  and 
depth-averaged  components  of  the  flow  and  salinity  as  dependent  variables. 


Conti nui t  y : 
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x-momentum : 
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y-mumentum : 
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Salinity : 
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An  equation  of  state  relating  the  water  density  to  the  salinity  and  water 
temperature  (assumed  constant)  can  be  written  as 


p(s,T)  =  1000  4-  ALO  *  P0 

where 

AL  =  1779.5  +  11.25T  -  0.0745T2  ~  (3.80  +  0.01T)s 
ALO  =  0.6980 

PO  =  5890.0  +  38T  -  0.375T2  +  3s 

In  the  above  equations,  the  surface  wind  shear  can  be  written  as 

w-  i 

c  2 

t  =  — -  p  v  cos  a 
s  p  aw 
X  o 


c  z 

t  =  —  p  v  sin  a 

s  paw 

y  o 


and  the  bottom  sh  ar  as 


l'he  Coriolis  parameter,  f  ,  is  computed  from 


f  =  2w  sin  X 
e 


(A16) 


where 

=  earth's  angular  velocity 

X  =  angle  of  latitude  of  the  center  of  the  area  being  modeled 
14.  In  order  to  finalize  the  above  system  of  equations,  it  remains  to 
couple  the  salinity  computations  with  those  of  the  flow  field.  Assuming  that 
the  pressure  is  hydrostatic,  it  can  be  shown  that  the  horizontal  pressure 
gradient  terms  can  be  written  as 


and 


Substituting 
final  form  of 
below: 
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Equations  A14-A17  into  Equations  A9-A11  would  then  yield  the 
the  governing  equations.  Variables  in  the  equations  are  defined 


=  Chezy  coefficient 

=  diagonal  components  of  eddy  viscosity  tensor 
=  off  diagonal  components  of  eddy  viscosity  tensor 
=  components  of  eddy  dispersion  tensor 
=  acceleration  due  to  gravity 
=  water  depth 
=  pressure 

=  atmospheric  pressure 
=  salinity 
=  temperature 
=  components  of  velocity 
=  wind  speed 

=  resistance  coefficient 
=  Cartesian  coordinates 


l  =  wind  direction 


>. 
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B 

x 


s 

X 


co  = 
e 

9/9 1  = 

3 _  3 _ 
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latitude  of  center  of  modeled  area 
water  density 
density  of  air 
reference  water  density 
components  of  bottom  shear  stress 

components  of  bottom  shear  stress 

water-surface  elevation 
earth's  angular  velocity 
time  derivative 

space  derivatives 


Two-Dimensional  Laterally  Averaged 


15.  The  equations  below  governing  flow  in  relatively  narrow  and  deep 
bodies  of  water  are  obtained  by  integrating  the  3-D  equations  over  the  width. 
Since  both  salinity  and  temperature  effects  may  influence  the  vertical  strati¬ 
fication,  a  transport-diffusion  equation  for  both  salt  and  temperature  is 
given. 


Longitudinal  (x-direction  momentum) 
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with 


t  =  C*p  /pW  cos  d>  (surface) 
z  a  a 
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( interlayer) 
(bottom) 
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Internal  continuity 


(wbb)k-l/2  =  (wbb)k+l/2  +  (UBh)  “  qBh 


Total  depth  continuity 


f-  -  X  ^ <uBh)  =IqBh 


Vertical  (z-direction)  momentum 
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Heat  balance 
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Salinity  balance 
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Equat ion  of  state 


p  =  [ 1000  P  /(LA  +  0.698  P  )] 
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where 


P  =  5890  +  38T  -  0.375T  +  3s 

o 


A9 


H  Bh 
n 

V 


sqBh 


LA  =  1779.5  +  11.25T  -  0.0745T  -  (3.8  +  0.01T)s 

Variables  in  the  equations  above  are  defined  below: 

2 

A  =  x-direction  momentum  dispersion  coefficient  (m  /s) 
x  2 

A^  =  z-direction  momentum  dispersion  coefficient  (m  /s) 
b  =  estuary  or  river  width  (m) 

B  =  laterally  averaged  width  integrated  over  h  (m) 

1/2 

C  =  Chezy  resistance  coefficient,  m  /s 
C*  =  resistance  coefficient  associated  with  wind 

2 

=  x-direction  temperature  and  salinity  dispersion  coefficient  (m  /s) 

2 

D  =  z-direction  temperature  and  salinity  dispersion  coefficient  (m  /s) 

z  2 
g  =  acceleration  due  to  gravity  (m/s  ) 

h  =  horizontal  layer  thickness  (m) 

=  source  strength  for  heat  balance  (C°m^s 

k  =  integer  layer  number,  positive  downward 

P  =  pressure 

3 

q  =  tributary  inflow  or  withdrawal  (m  /s) 
s  =  laterally  averaged  salinity  integrated  over  h  (ppt) 
t  =  time  (s) 

T  =  laterally  averaged  temperature  integrated  over  h  (C°) 

U  =  x-direction,  laterally  averaged  velocity  integrated  over  h  (m/s) 

u  =  x-direction,  laterally  averaged  velocity  (m/s) 

3 

V  =  cell  volume  (B  •  h  •  Ax)  (m  ) 

W  =  wind  speed  (m/s) 

3 

w^  =  z-direction,  laterally  averaged  velocity  (m/s) 

x,  z  =  Cartesian  coordinates:  x  is  along  the  estuary  center  line  at 
the  water  surface,  positive  to  the  right  and  z  is  positive 
downward  from  the  x-axis  (m) 

>,  =  surface  elevation  (m) 

3 

p  =  density  (kg/m  ) 

P  =  air  density  (kg/m^) 

3 

i  =  z-direction  shear  stress 
z 

}>  =  wind  direction  (rad) 


16.  The  concept  of  eddy  coefficients  is  utilized  to  represent  the  effect 
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of  both  time  averaging  and  spatial  averaging  of  the  equations.  The  horizontal 
coefficients,  and  ,  are  assumed  to  be  constant;  whereas  the  vertical 

coefficients,  and  ,  are  dependent  upon  the  stratification  as  re- 

"lectcd  by  the  local  Richardson  number,  ,  i.e., 

A  =  A  / 1  +  3.33R. 

Z  Z  \  X 

O  '■ 


where 


(A25) 


(A26) 


and  A  and  D  are  the  vertical  coefficients  for  no  stratification.  Due 

z  z 

to  the  hydrostatic  pressure  assumption,  unstable  stratification  cannot  be 
modeled  in  a  convective  fashion  and  thus  is  handled  in  a  diffusive  manner  by 
increasing  to  its  stability  limit. 

17.  The  laterally  averaged  horizontal  pressure  gradient  in  the  longi¬ 
tudinal  momentum  equation  contains  the  density  driving  force.  Using  the 
expression  for  the  hydrostatic  pressure,  the  horizontal  pressure  gradient  can 
be  divided  into  its  two  components  of  the  barotropic  (surface  slope)  gradient 
and  the  baroclinic  (density)  gradient  to  yield: 


(A27) 


Substituting  liquations  A25-A27  into  Equations  A18  and  A22-A23  would  then  yield 
the  final  form  of  the  governing  equations  for  the  computation  of  laterally 
averaged  stratified  flows. 
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